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A  bstract 

Accurate  numerical  modeling  of  complex,  multi-dimensional  shock 
propagation  is  needed  for  many  Department  of  Defense  applications.  A  three- 
dimensional  code,  based  upon  E.F.  Toro’s  weighted  average  flux  (WAF)  method 
has  been  developed,  tested,  and  validated.  Code  development  begins  with  the 
introduction  and  application  of  all  techniques  in  a  single  dimension.  First-order 
accuracy  is  achieved  via  Godunov’s  scheme  using  an  exact  Riemann  solver. 
Adaptive  techniques,  which  employ  approximate  solutions,  are  implemented  to 
improve  computational  efficiency.  The  WAF  method  produces  second-order 
accurate  solutions,  but  introduces  spurious  oscillations  near  shocks  and  contact 
discontinuities.  Total  variation  diminishing  (TVD)  flux  and  weight  limiting 
schemes  are  added  to  reduce  fluctuation  severity.  Finally,  the  fully  developed 
one-dimensional  code  is  validated  against  experimental  data,  and  extended  into 
two  and  three  dimensions  via  dimension-splitting  techniques. 
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DEVELOPMENT  OF  A  THREE-DIMENSIONAL  AIR  BLAST 
PROPAGATION  MODEL  BASED  UPON  THE  WEIGHTED  AVERAGE 

FLUX  METHOD 

1.  Introduction 

1.1.  Motivation 

The  shock  wave  is  an  important  natural  phenomenon  which  commands 
attention.  Accordingly,  accurate  and  efficient  numerical  modeling  of 
compressible,  inviscid,  time-dependent  flow  is  a  major  branch  of  computational 
fluid  dynamics.  Analysis  of  complex,  multi-dimensional  shock  propagation  is  an 
application  of  this  numerical  modeling  and  is  needed  for  many  Department  of 
Defense  (DoD)  applications. 

1.2.  Background 

Fluid  dynamics  is  governed  by  a  system  of  non-linear,  hyperbolic  partial 
differential  equations  which  dictate  conservation  laws  of  mass,  momentum,  and 
energy.  In  1952,  Courant,  Isaacson,  and  Rees  introduced  a  first-order  accurate 
scheme  to  solve  these  equations  by  finite  differences.  In  1959,  Sergei  K.  Godunov 
presented  an  extension  to  this  scheme  which  resolved  a  computational  domain  as 
a  piecewise  constant  distribution  of  the  conserved  variables  in  a  series  of 
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boundary  value  problems  (Godunov,  1959:1187).  These  boundary  value 
problems,  known  as  Riemann  problems,  are  solved  exactly  to  satisfy  the 
governing  conservation  laws  and  determine  numerical  fluxes  of  mass,  momentum, 
and  energy  at  each  boundary.  Godunov’s  method  then  applies  these  boundary 
fluxes  to  calculate  the  conserved  variable  values  at  either  side  of  the  boundary, 
which  serve  as  conditions  to  solve  Riemann  problems,  forward  in  time,  as 
required. 

Due  to  the  non-linearity  of  the  conservation  laws,  numerical  iteration  is 
required  to  obtain  exact  solutions  to  the  Riemann  problems.  Since  the 
introduction  of  Godunov’s  method,  several  approximate  Riemann  solvers  have 
been  developed  that  predict  solutions  without  requiring  computationally  costly 
iteration.  Unfortunately,  Godunov’s  method,  implemented  with  exact  or 
approximate  Riemann  solvers,  provides  only  first-order  accuracy,  making  it 
unsuitable  for  application  to  practical  problems  where  well-resolved  solutions 
require  fine  meshes.  Specifically,  first-order  schemes  often  result  in  smeared 
solutions  at  discontinuities,  such  as  shock  fronts,  and  predict  slow  convergence  to 
these  solutions  (Toro,  1999:213). 

In  his  text,  Riemann  Solvers  and  Numerical  Methods  for  Fluid 
Dynamics,  2d  Edition ,  E.F.  Toro  presents  high-resolution,  second-order 
techniques  to  remedy  the  shortcomings  associated  with  first-order  methods. 
Specifically,  Toro’s  techniques  permit  determination  of  more  accurate  numerical 
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fluxes  at  each  boundary  than  Godunov’s  method  allows.  These  fluxes  are 
calculated  by  analyzing  the  wave  structure  of  the  Riemann  problem  to  determine 
a  weighted  average  of  the  flux  vectors  at  the  boundary.  While  Toro’s  weighted 
average  flux  (WAF)  method  improves  accuracy,  it  also  introduces  false 
fluctuations  at  discontinuities,  such  as  shock  fronts.  Toro  ultimately  applies  total 
variation  diminishing  (TVD)  techniques  to  limit  the  flux  weights  and  resolve 
these  discontinuities  in  the  regions  of  interest.  Finally,  Toro’s  methods  may  be 
applied  to  solve  single  or  multi-dimensional  problems  (Toro,  1999:490). 

1.3.  Problem 

The  principal  endeavor  of  this  research  is  the  development  of  a  three- 
dimensional  hydrodynamic  shock  code  to  solve  the  Euler  conservation  equations 
based  upon  the  weighted  average  flux  method  (Toro,  1999:492).  The  secondary 
objective  applies  the  developed  code  to  model  air  blast  propagation  in  one,  two, 
and  three  dimensions. 

1.4 •  Scope 

The  developed  computational  model  is  limited  to  modeling  air  blast  under 
ideal  conditions. 

1.5.  Assumptions  and  Limitations 

The  numerical  hydrodynamics  code  assumes  air  behaves  as  an  ideal  gas 
with  a  constant  ratio  of  specific  heats  (y=  1.4).  This  assumption  limits 
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application  of  the  code  to  low  and  intermediate  pressure  regions  where  peak 
pressures  remain  below  100  psi  (690  kPa).  Modeling  of  air  blast  propagation 
assumes  that  the  atmosphere  is  homogeneous,  with  sea  level  properties,  up  to 
5000  feet  (Glasstone  and  Dolan:  103).  The  model  also  assumes  the  earth  behaves 
as  an  ideal  surface  which  is  flat  and  reflects  all  incident  energy. 

1.6.  Approach 

Accurate  and  efficient  one-dimensional  numerical  hydrodynamics  modeling 
serves  as  the  basis  for  code  development  in  two  and  three  dimensions.  An 
algorithm  to  solve  a  single  Riemann  problem  for  the  Euler  conservation 
equations,  exactly,  was  developed  first.  Next,  a  Godunov’s  method  algorithm 
applied  the  exact  Riemann  solver  to  an  entire  computational  domain,  providing 
first-order  accuracy.  Godunov’s  algorithm  then  incorporated  three  approximate 
Riemann  solvers,  in  order  to  maximize  computational  efficiency.  Application  of 
Toro’s  WAF  method  followed  and  provided  second-order  accurate  solutions, 
while  four  TVD  weight  limiters  were  incorporated  to  reduce  false  fluctuations  at 
discontinuities.  During  each  step  of  code  development,  numerical  values  of  the 
physical  variables  of  density,  material  speed,  and  pressure,  as  well  as  internal 
energy,  were  verified  against  five  documented  shock  tube  tests  using  a  fine  mesh 
of  10001  points  over  a  domain  length  of  one  meter  for  the  exact  Riemann  solver, 
and  a  coarse  mesh  of  200  points  for  each  of  the  Godunov  method  algorithms. 
Finally,  the  second-order  accurate  Godunov  WAF  TVD  model  was  validated 
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against  the  Army  Research  Laboratory  (ARL)  57cm  shock  tube  test  using  a  mesh 
of  5000  points  over  50  meters. 

Following  validation  of  the  one-dimensional  model,  the  code  was  extended 
into  two  and  three  dimensions  by  dimension  splitting  techniques.  The  two  and 
three-dimensional  codes  were  verified  against  mild  cylindrical  and  spherical  shock 
tests  presented  by  Toro. 
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2.  Fluid  Dynamics  Theory 


2.1.  Governing  Equations 

The  laws  of  conservation  of  mass,  momentum,  and  energy  as  applied  to 
fluid  motion  are  summarized,  comprehensively,  by  a  set  of  non-linear  partial 
differential  equations  known  as  the  Navier-Stokes  equations.  For  mathematical 
and  experimental  reasons,  the  boundary  conditions  which  accompany  these 
computationally  precise  equations  dictate  solution  procedures  which  are  often 
inefficient  and  impractical.  However,  application  of  some  simplifying 
assumptions  reduces  these  complex  equations  to  a  simpler  set  of  partial 
differential  equations,  commonly  referred  to  as  the  Euler  equations  (Chorin  and 
Marsden,  1990:34). 

2.2.  Simplifying  Assumptions 

The  Euler  equations  result  from  the  assumptions  that  viscosity,  heat 
conduction,  and  body  forces  are  negligible  for  a  compressible  medium.  If  these 
assumptions  are  maintained  by  the  computational  model  that  simulates  shock 
propagation  through  air  and  shock  reflection  at  a  surface,  then  the  Euler 
equations  may  be  justifiably  applied  to  this  model.  Therefore,  an  analysis  of 
viscosity,  heat  conduction,  and  body  forces  within  the  air,  at  the  boundary 
between  the  air  and  surface,  and  at  the  surface  where  shock  interaction  occurs  is 
necessary. 
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The  computational  model  assumes  a  homogeneous  atmosphere  with  sea- 


level  properties.  Regarding  air  located  far  from  a  boundary,  this  assumption 
implies  negligible  density,  velocity,  and  temperature  gradients  exist,  and  justifies 
ignoring  body  forces,  viscosity,  and  heat  conduction,  respectively. 

The  model  also  assumes  that  the  earth  behaves  as  an  ideal  surface  that  is 
flat  and  reflects  all  incident  energy,  which  allows  heat  conduction  between  the  air 
and  the  surface  to  be  ignored.  Furthermore,  shocks  propagating  over  this 
smooth,  flat  surface  exert  high  pressures  and  inertias  over  a  physical  domain  that 
is  large  in  comparison  to  the  thin  surface  boundary.  Therefore,  variations  in 
pressure  normal  to  the  surface  are  negligible  within  the  boundary  layer,  which 
indicates  the  pressure  distribution  at  the  surface  results  from  the  flow  above  the 
boundary.  Hence,  viscosity  at  the  boundary  between  the  air  and  surface  can  also 
be  ignored.  Finally,  because  viscosity  and  heat  conduction  may  be  ignored 
within  air  located  far  from  a  boundary,  and  at  the  boundary  between  the  air  and 
the  ideal  surface,  the  air  at  these  locations  is  characterized  as  inviscid  and 
adiabatic  (Wittig,  1999:5). 

Examination  of  air  located  at  the  shock  front  indicates  this  air  is  neither 
inviscid  nor  adiabatic.  Physically,  a  shock  front  represents  a  thin  region  where 
the  properties  of  density,  pressure,  and  material  velocity  change  rapidly.  The 
temperature  and  velocity  gradients  incurred  at  the  front  imply  viscosity,  heat 
transfer,  and  entropy  cannot  be  ignored,  and  the  simplifying  assumptions  that 
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allow  reduction  of  the  Navier- Stokes  equations  to  the  Euler  equations  are 
effectively  threatened.  However,  if  the  physical  shock  front  is  represented, 
mathematically,  as  a  discontinuity,  and  the  width  of  this  discontinuity  is 
assumed  to  be  small  with  respect  to  the  domain  represented  by  the 
computational  model,  and  the  second  law  of  thermodynamics  sufficiently 
recognizes  the  discontinuity  and  predicts  physically  correct  behavior  at  the  shock 
front,  then  viscosity  and  heat  conduction  may  be  neglected  at  the  front,  as  an 
approximation.  In  fact,  smooth  solutions  obtained  from  the  Navier-Stokes 
equations  compared  to  discontinuous  solutions  obtained  by  applying  simplifying 
assumptions  are  nearly  indistinguishable  (LeVeque,  1992:9). 

The  assumptions  that  air  viscosity,  heat  conduction,  and  body  forces  are 
negligible  justify  applying  the  Euler  equations  to  the  developed  computational 
model. 

2.3.  The  Euler  Equations 

The  time-dependent  Euler  form  of  the  conservation  equations  of  mass, 
momentum,  and  energy  may  be  written  in  compact  vector  notation  as 

U,  +  FtU)*  +  G(U)„  +  H(U),  =  0  (2.1) 

where  U  is  a  column  vector  of  conserved  variable  densities  and  F(U),  G(U),  and 
H(U)  represent  fluxes  of  conserved  variables  in  the  x,  y,  and  z  directions, 
respectively,  as 
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Vectors  are  denoted  by  boldface  letters.  The  subscript  t  in  equation  (2.1) 
represents  a  temporal  partial  derivative,  while  subscripts  x,  y,  and  2  symbolize 
spatial  partial  derivatives  (Toro,  1999:3). 

Two  sets  of  variables,  conserved  and  primitive,  are  used  in  the 
conservation  equations,  and  there  exists  some  freedom  when  choosing  variables  to 
describe  the  flow  under  consideration.  Conserved  variables  represent  conserved 
quantities  of  mass,  momentum,  and  energy,  while  the  primitive  variables  of 
density,  velocity,  and  pressure  can  be  easily  measured  and  are  often  imposed  by 
the  domain  of  the  problem.  Algebraic  simplification  and  application  of  the  ideal 
gas  equation  of  state,  discussed  in  Section  2.4,  allows  the  primitive  variables  to 
be  expressed  in  terms  of  the  conserved  variables,  or  vice  versa,  where 
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In  a  single  dimension,  assuming  uniform  flows  in  the  y  and  z  directions, 
equation  (2.1)  may  be  written 

Ut  +  F(U)I  =  0  (2.5) 

where 


p 

pu 

u  = 

pu 

,  F  = 

to 

+ 

E 

u(E  +  p) 

(2.6) 


Similarly,  a  two-dimensional  analysis  where  uniform  flow  exists  in  the  .2-direction 
requires  solving 

U,  +  FfU),  +  G(U)„  =  0  (2.7) 

where 


U 


P 

pu 

pv 

pu 

to 

+ 

pvu 

pv 

,  F  = 

,  G  = 

2  , 

puv 

pv  +  p 

E 

u(E  +  p) 

v{E  +  p) 

(2.8) 
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Examination  of  the  single  or  multi-dimensional  Euler  equations  reveals,  in  all 


cases,  that  the  number  of  unknown  variables  exceeds  the  number  of  equations  by 
one.  Therefore,  an  appropriate  thermodynamic  equation  of  state  relating  two  or 
more  physical  quantities  within  the  medium  is  required  to  guarantee  a  solution 
for  all  unknown  variables. 

2-4-  Thermodynamic  Considerations 

The  energy  density,  E,  included  in  equation  (2.2)  is  expressed 

E  =  —p(u2  +  v2  +  w2  )  +  pe  (2.9) 

2 

where  e  is  specific  internal  energy.  Assuming  the  medium  within  the 
computational  model  behaves  as  a  calorically  ideal  gas,  then  the  density  and 
pressure  within  the  medium  are  related  by  the  ideal  gas  equation  of  state  (EOS): 

p  =  pRT  (2.10) 

where  R  is  the  universal  gas  constant  and  T  is  the  temperature.  If  the  ideal  gas 
EOS  is  assumed,  then  it  follows  that  internal  energy  is  a  function  of  temperature 
alone,  or 

e  =  cvT  (2.11) 

where  cv  ,the  specific  heat  capacity  at  constant  volume,  is  a  constant. 
Furthermore, 
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cv  —  °v  —  R 


(2.12) 


where  cp  ,the  specific  heat  capacity  at  constant  pressure,  is  also  constant. 
Defining  the  ratio  of  specific  heats,/,  as 


(2.13) 


and  substituting  equations  (2.11)  and  (2.12)  into  equation  (2.10)  yields  the 


relationship 

e  = - — - . 


(2.14) 


Substitution  of  (2.14)  into  (2.9)  presents  a  relationship  for  E  in  terms  of  the 
primitive  variables  p,  u,  v,  w,  and  p  and  provides  necessary  closure  to  the  system  of 

Euler  conservation  equations  (Toro,  1999:13). 
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3.  Exact  Riemann  Solver  Development 


3.1.  The  Riemann  Problem  for  the  Euler  Equations 


The  Riemann  problem  (RP)  for  the  one-dimensional  time-dependent  Euler 


conservation  equations  is  the  initial  value  problem  (IVP) 


U,  +  F(U),  =  0 

p 

pu 

u  = 

pu 

,  F  = 

£ 

to 

+ 

5 

E 

u(E  +  p ) 

(3.1) 


with  initial  conditions 


U(x,  0)  =  U{a\x) 


UL 

x  <  0 

UR 

x  >  0 

(3.2) 


The  domain  of  interest  in  the  x-t  plane  is  the  set  of  points  (x,  t )  where 
-oo  <  x  <  oo  and  t>  0.  In  practice,  x  varies  in  a  finite  interval  [xL,xR\  around  the 
point  x  =  0.  Additionally,  the  vector  of  primitive  variables  at  a  specified  time 
W  =  ( p,u,p )  is  frequently  used  when  solving  the  RP,  as  opposed  to  the  vector  U 

of  conserved  variables.  Therefore,  the  RP  is  the  IVP  for  (3.1)  consisting  of 
constant  data  states  WLandWR  left  and  right  of  a  discontinuity  located  at  x  = 

0. 

Physically,  as  applied  to  the  Euler  equations,  the  RP  represents  a 
simplified  shock  tube  problem.  In  this  problem,  two  stationary  gases  are  initially 
separated  by  a  boundary.  Removal  of  this  boundary  generates  a  nearly  centered 
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system  of  three  waves  that  typically  consists  of  a  rarefaction  wave,  contact 
discontinuity,  and  a  shock  wave  (Toro,  1999:117;  Ben-Artzi,  2003:98). 

3.2.  The  Riemann  Problem  Wave  Structure 

Figure  3.1  depicts  one  example  of  a  Riemann  problem  wave  structure. 
Solutions  to  the  system  of  partial  differential  equations  in  the  IVP  (3.1)  produce 
three  wave  families  associated  with  the  eigenvalues  listed  above  each  wave 
number  in  Figure  3.1.  For  each  eigenvalue,  u  is  the  material  velocity  and  a 
represents  the  local  speed  of  sound. 


Figure  3.1.  Riemann  Problem  Wave  Structure 

As  indicated  in  the  figure,  the  wave  families  effectively  separate  the  RP,  from  left 
to  right,  into  four  regions  Wl,  W*l,  W*jj,,  and  Wp  .  Although  the  left  and  right 
most  states  are  known  based  on  the  initial  conditions  specified  in  the  RP, 
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information  about  the  star  region  is  necessary  in  order  to  fully  define  the  RP. 
Solutions  dictate  waves  one  and  three  may  be  classified  as  shocks  or  rarefactions, 
while  wave  two  is  always  a  contact  surface  wave.  Traditionally,  a  rarefaction  is 
represented  as  a  pair  of  rays,  with  a  head  and  tail,  while  shocks  are  depicted  by  a 
single,  solid  line.  Until  the  full  classification  of  the  outer  waves  as  shocks  or 
rarefactions  occurs,  both  waves  are  portrayed  as  rarefactions.  Consequently, 
accurate  identification  of  waves  one  and  three  as  shocks  or  rarefactions  must 
occur  before  full  definition  of  the  RP,  and  the  method  of  characteristics  provides 
the  ability  to  accomplish  this  identification. 

3.3.  The  Method  of  Characteristics 

Characteristic  vectors,  often  called  eigenvectors,  provide  straight-line 
solutions  to  a  system  of  linear  differential  equations.  Analysis  of  these  vectors 
facilitates  generation  of  characteristic  fields,  which  provide  qualitative 
information  about  the  behavior  of  a  system  of  equations  (Blanchard,  2002:258). 
Unfortunately,  the  Euler  equations  are  non-linear,  resulting  in  discontinuities 
where  characteristics  cross  (Wittig,  1999:14).  Nevertheless,  the  method  of 
characteristics  generates  useful  qualitative  information  and  accurately  identifies 
RP  waves  as  shocks,  rarefactions,  or  contact  surface  waves.  This  classification, 
coupled  with  known  information  from  the  left  and  right  states  specified  by  the 
initial  conditions  facilitates  solution  for  variables  in  the  star  region  and, 
ultimately,  completes  solution  for  variables  in  all  four  RP  regions. 
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Recall  that  properties  of  density,  material  velocity,  and  pressure  are 


mathematically  discontinuous  for  a  shock  wave.  Conversely,  a  rarefaction  wave 
is  a  smooth  wave  region  where  solutions  for  the  same  properties  vary 
continuously  between  the  head  and  tail.  A  characteristic  vector  field  associated 
with  a  right  moving  rarefaction  is  depicted  in  Figure  3.2.  The  divergent 
characteristics  within  the  fan  and  on  either  side  of  the  wave  indicate  pressure  and 
density  decrease  across  the  expansion  fan  from  head  to  tail,  while  material 
velocity  also  decreases.  A  left  moving  rarefaction  exhibits  identical  behavior, 
except  for  material  velocity,  which  increases  across  the  expansion  fan  from  head 
to  tail.  In  both  cases,  the  speed  of  the  head  represents  the  speed  of  the 
rarefaction  wave. 

t 


Figure  3.2.  Characteristic  Field  of  a  Right  Moving  Rarefaction 

Figure  3.3  depicts  a  characteristic  vector  field  associated  with  a  left 
moving  shock.  Along  the  line  where  characteristics  cross,  solutions  are 
discontinuous.  Accordingly,  properties  of  density,  material  velocity,  and  pressure 
experience  an  instantaneous  jump  across  the  wave  from  head  to  tail. 
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Graphically,  the  head  and  tail  are  conventionally  removed  and  the  shock  wave  is 
represented  by  a  bold,  solid  line.  Finally,  the  speed  of  the  head  does  not 
represent  shock  speed,  as  the  shock  wave  cannot  be  classified  characteristic.  This 
shock  speed  will  be  analyzed  later. 


t 


Figure  3.3.  Characteristic  Field  of  a  Left  Moving  Shock 

Recall  the  center  wave  in  the  RP  structure  is  always  a  contact  surface 
wave.  Figure  3.4  shows  the  characteristic  vector  field  is  parallel  to  the  contact 
surface. 


Figure  3.4.  Characteristic  Field  of  a  Right  Moving  Contact  Surface 

In  the  context  of  the  RP,  both  pressure  and  material  velocity  are  constant  across 
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wave  two  as 

P*  L  =  P*  R  =  P*> 

U*i^  =  U*  =  u*  7 

while  density  experiences  a  jump  across  the  contact  surface  and 

P*L^P*R-  (3-4) 

Classification  of  waves  one  and  three  as  shocks  or  rarefactions,  together 
with  primitive  variable  information  from  states  WL  and  WR  provides  necessary 
information  to  solve  forW*L  and  W*R  and  fully  define  the  RP.  The  exact 
Riemann  solver  (ERS)  provides  the  method  to  accomplish  this  task. 

3.4 -  The  Exact  Riemann  Solver  (ERS) 

3-4-1 -  Solution  Strategy 

As  a  pedagogical  exercise,  Toro  applies  the  ERS  to  a  simplified  shock  tube 
problem.  Specifically,  the  shock  tube  is  treated  as  a  single  RP,  with  a  diaphragm 
separating  constants  states  left  and  right  of  the  boundary.  Methods  for 
determining  the  values  of  primitive  variables  in  the  star  region  are  presented. 
Following  solution  for  these  variables  in  all  four  RP  states,  primitive  variable 
values  are  calculated  at  a  user  specified  number  of  points  within  the  shock  tube. 
Toro  refers  to  this  second  step  as  sampling.  While  this  method  is  a  simplified 
approach  to  shock  tube  modeling,  the  solution  of  a  single  RP  serves  as  the  basis 


3.3) 
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for  subsequent  methods  which  treat  a  computational  domain  as  a  series  of  local 


Riemann  problems. 

3-4-2.  Calculation  of  Star  Variables 

Toro’s  process  to  connect  known  and  unknown  states  is  straight-forward. 
Material  velocities  in  the  left  and  right  star  states  are  determined  using  the 
equations 

U*L  =  UL  -fL(p*,WL), 

(3. 

U*R  =  UR  +  fR(p*  ,WL)  . 


The  functions  /L  and  /R  are  dependent  on  the  classification  of  the  left  and 
right  waves.  If  the  outer  wave  is  a  shock,  the  Rankine-Hugoniot  conditions  are 
applied  to  determine  the  functions.  On  the  other  hand,  if  the  outer  wave  is  a 
rarefaction,  isentropic  relationships  and  the  Generalized  Riemann  Invariants 
along  the  characteristics  are  used  (Toro,  1999:120).  A  generalized  pressure 
function, W^) ,  where  E,  indicates  the  appropriate  left  or  right  state,  results 
as 


/{(?*,  w4) 


(p*  -  !>|) 


2  at 


7-1 


(  Y 


+  -0  +  P?(7  -1)) 

7-1 


if  p*  >  (shock) , 


p* 


27 


1 


(3,6) 


if  p*  <p^  (rarefaction). 
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It  follows  from  (3.3)  that  within  the  star  region  u *L  —  u*R  =  0 .  Applying  this 
relationship  to  (3.5)  yields  the  equation 

fl(P*>  Wl)  +  Ir(p*,  Wr)  +  ur  —ul  =0  (3.7) 

whose  root,  p* ,  must  be  found.  The  process  for  determining  numerically  p*  is,  at 
this  point,  uncomplicated.  First,  an  initial  guess  for p* is  generated.  Next,  the 
generated  value  is  substituted  into  (3.6)  in  order  to  determine  fL(p*,WL) 
and  fR(p* ,W R) .  Finally,  the  condition  defined  in  (3.7)  is  verified  against  a 
specified  tolerance.  If  the  condition  is  satisfied  the  process  is  complete,  but  if  the 
condition  fails,  the  process  is  repeated. 

Newton’s  method  was  chosen  as  the  root  finding  algorithm  due  to  its 
quadratic  convergence  to  a  solution  and  ease  of  implementation.  Toro  facilitates 
incorporating  this  algorithm  by  providing  explicit  derivatives  for  the  pressure 
functions,  which  Newton’s  method  requires,  and  by  presenting  techniques  which 
generate  initial  guesses  for  p* ,  based  on  known  primitive  variable  information 
from  the  left  and  right  states. 

After  determining  pressure  in  the  star  region,  the  solution  ioru*  follows 
(Toro,  1999:119)  as 

u*  =  ~(uL  +  uR )  +  ~[fR  (p*,WR)  +  fL  (p*,WL)\  (3.8) 

Recall  within  the  star  region  p*L  ^  p*R .  The  values  for  density  within 
this  region  depend  on  the  classification  of  the  wave  and  are  determined  via  the 
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generalized  density  function  pk^(p*,p^,p^) ,  where  %  indicates  the  appropriate  left 
or  right  state  (Toro,  1999:133),  as 


P*s(p*,P£,At) 


Pi 


P*(7  +  l)  +  ^(7~l) 

p*h~1)  +  Pih  +  1))/ 


if  p*  >  p ^  (shock) , 


Pi 


if  P*  <Pi  (rarefaction). 


(3.9) 


Solution  for  primitive  variables  in  the  star  region  is  complete,  and  the  RP  is  fully 
defined. 

3.f..3.  S  amp  ling 

Sampling  determines  the  exact  solution  to  the  RP  at  any  point  (x,  t )  in 
the  domain  of  interest.  The  procedure  makes  use  of  an  important  mathematical 
property  that  states,  along  straight  rays  emanating  from  the  origin  (x  =  0), 
solutions  to  the  RP  are  self-similar  (Ben-Artzi,  2003:99).  Solutions  at  different 
time  instants  are  obtained  from  one  another  by  a  similarity  transformation, 
which  is  analogous  to  scaling.  In  short,  the  knowledge  of  the  RP  solution  at  an 
instant  in  time  t0  is  sufficient  to  obtain  a  solution  for  all  t  >  0  (Polyanin, 
2004:695). 

The  physical  domain  for  the  shock  tube  problem  is  [0,  L]  x  [0,  T]  where  T 
is  the  specified  solution  time.  To  guarantee  self-similarity  within  the  problem, 
the  discontinuity,  physically  located  at  x  =  L/ 2  in  the  shock  tube  coordinate 
system,  must  be  transformed  to  x'  =  0 ,  where  x'  represents  the  location  in  the 
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RP  coordinate  system.  Figure  3.5  illustrates  that  the  coordinate  transformation 


occurs  for  all  points  x'  as 


x 


L 

x - . 

2 


(3.10) 


It  follows  that  the  speed  required  to  reach  a  point  x1  at  time  t  in  the  RP 
coordinate  system  is 

,'  =  ^.  (3.11) 


Figure  3.5a.  Shock  Tube  Coordinates  Figure  3.5b.  RP  Coordinates 


A  solution  for  the  primitive  variables  at  the  point  (x1,  t)  is  obtained  by 
comparing  s'  to  the  three  RP  wave  speeds.  These  wave  speeds  are  determined 
by  the  set  of  equations 
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1 

2 


if  p*  >  pL  (shock), 
if  p*  <  pL  (rarefaction), 


7  +  1 

/  \ 

P* 

Si=- 

uL  —  aL 

27 

.Pl, 

+  7^1 
27 


UL  ~  aL 


S2  —  U*, 


(3.12) 


^3 


UR  +  aR 


7  +  1 

/  \ 

P* 

7-1 

27 

kPR  j 

27  . 

UR  +  aR 


if  p*  >  pR  (shock), 
if  p*  <  pR  (rarefaction). 


If  a  left  rarefaction  is  present,  then  the  speeds  of  the  head  and  tail  are  calculated 
as 


Shl  —  uL  —  «l  • 

StL  =  U*L  —  a*L  ■ 

For  a  right  rarefaction,  the  head  and  tail  move  with  speeds 


(3.13) 


Shr  —  UR  +  aR, 


StR  —  U*R  +  a*R> 


(3.14) 


where 


a*t 


ip* 


{ H 


(3.15) 


At  this  point,  all  necessary  information  to  determine  W(x',  t )  exists. 
Because  u*  represents  the  contact  surface  wave  speed,  comparison  of  s'  to  u* 
determines  if  the  sample  point  is  located  left  or  right  of  the  contact  discontinuity. 
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Specifically,  if  s'  <  u*  the  sample  point  is  located  left  of  the  contact. 

Conversely,  the  point  is  right  of  the  contact  wave  when  s'  >u* . 

Figure  3.6  shows  when  the  sample  point  (x',  t )  lies  left  of  a  contact  discontinuity, 
and  in  the  presence  of  a  shock,  primitive  variables  are  determined  as 

WL  if  —  <  SL, 

W  (x\t)  =  -  1  ,  (3.16) 

W^L  if  SL  <  j  <  u* . 

On  the  other  hand,  in  the  presence  of  a  rarefaction,  primitive  variables  are 
calculated  as 


Figure  3.6.  Sampling  a  Solution  Point  Left  of  the  Contact 
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The  same  analysis  occurs  when  determining  primitive  variable  values  for  a 
point  (x1,  t )  right  of  a  contact  surface  discontinuity.  Figure  3.7  shows,  in  the 
presence  of  a  shock,  primitive  variables  are  calculated  by 


Figure  3.7.  Sampling  a  Solution  Point  Right  of  the  Contact 
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In  the  presence  of  a  rarefaction,  primitive  variables  are  determined  as 

W*R  if  u*<  —  < 

W(V,  t)  =  •  W Rfan  if  STR  <  —  <  5hk  ,  (3.20) 

Wr  if  —  >  6'hr  ; 

where  variables  within  the  fan  are  calculated  (Toro,  1999:136) 

2 

7-1 

; 

(3. 
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7-1 

The  entire  solution  of  the  Riemann  problem  is  defined. 

3. 5.  Validation  of  the  ERS  -  The  Sod  Shock  Tube  Test 

The  developed  Riemann  solver  was  validated  against  a  classical  shock  tube 
problem  presented  by  Sod  in  1978.  The  shock  tube  is  one  meter  long,  and 
contains  a  diaphragm,  positioned  at  x  =  0.5m,  which  separates  two  constant 
initial  states 

WL  =  (1.000  kg/m3, 0.0  m/sec,  1.0  N/m2)  x  <  0.5m, 

W(x,  0)  =  \ 

Wr  =  (0.125  kg/m3, 0.0  m/sec,  0.1  N/m2)  x  >0.5m. 

The  solution,  consisting  of  a  left  rarefaction,  contact,  and  right  shock  wave 
pattern,  develops  over  the  time  interval  t  =  [0.0,  0.25]  sec.  Profiles  of  primitive 


2  7  —  1  ,  x' 

P  PR7  +  1  aR(7  +  l)  t 


WR  fanO^R,  x'i  t)  ~  \  u  ~ 


2  [  (7  —  1)  X] 

U  — - (In  H - Un  H - 

7  + 1  [  2  K  t 


2  7  —  1  ,  x\ 
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variables  and  internal  energy  at  10,001  equally-spaced  points  throughout  the 


shock  tube  are  depicted  in  Figure  3.8  -  Figure  3.11.  In  addition  to  the  Sod  test, 
four  other  shock  tests  were  solved  using  the  ERS.  The  solutions  to  these  tests 
consisted  of  combinations  of  strong  shocks  and  rarefactions  (Einfeldt  et.  al., 
1991:273-295;  Woodward  and  Colella,  1984:115-173).  Numerical  results  proved 
consistent  with  published  data  for  all  tests. 

The  exact  Riemann  solver  algorithm  serves  as  the  basis  for  subsequent 
methods  which  treat  a  computational  domain  as  a  series  of  local  Riemann 
problems.  Furthermore,  the  solver  applied  to  the  Sod  shock  tube  test  provides  a 
benchmark  for  comparing  calculations  determined  by  finite  volume  techniques. 
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Figure  3.8.  Density  Profile  at  t  =  0.25  sec 
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Figure  3.9.  Material  Velocity  Profile  at  t  =  0.25  sec 
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Pressure  (N/m2) 
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Figure  3.10.  Pressure  Profile  at  t  =  0.25  sec 
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4-  Code  Implementation  in  One  Dimension 


4.1.  Overview 

A  general  scheme  for  non-linear  hyperbolic  systems  of  equations  is 
presented  in  this  chapter.  Finite  volume  techniques  are  introduced  and  applied 
to  the  time-dependent  Euler  conservation  equations  in  a  single  dimension. 
Godunov’s  first-order  upwind  conservative  method,  which  employs  exact 
solutions  to  local  Riemann  problems  in  order  to  calculate  numerical  intercell 
fluxes,  is  described.  Improvements  to  computational  efficiency  are  introduced  as 
approximate  Riemann  solvers.  Second-order  extensions  of  Godunov’s  method 
and  total  variation  diminishing  (TVD)  techniques  are  implemented  to  improve 
solution  accuracy.  The  developments  presented  allow  relaxation  of  symmetry, 
which  Toro’s  ERS  assumes,  in  order  to  obtain  solutions  to  more  extensive 
categories  of  problems  than  the  ERS  permits.  The  developed  computational 
model  is  validated  against  the  benchmark  shock  tube  test,  as  well  as  a 
documented  experimental  shock  test. 

4-2.  Initial  Boundary  Value  Problem  (IBVP) 

The  one-dimensional  initial  boundary  value  problem  for  non-linear  systems 
of  hyperbolic  conservation  laws,  assuming  uniform  flows  in  the  y  and  z  directions 
is  defined  as 
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PDEs  :  Vt  +  F(U)a;  =  0, 

ICs:  U(x,0)  =  U(0)(a:), 

BCs  :  U(0,  t)  =UL(0;  U (L,  t )  =  U R(t). 


(4.1) 


U(x,  t)  is  a  vector  of  conserved  variables;  F(U)  is  a  vector  of  fluxes;  U1(^(x)  is  the 
piecewise  constant  distribution  of  initial  data  at  time  t  =  0;  [0,  L]  is  the  spatial 
domain  and  the  boundary  conditions  are  represented  by  boundary  functions 
UL(f)and  UR(t)(Toro,  1999:213).  Boundary  conditions  may  be  specified  as 
symmetric,  or  reflective,  or  transmissive. 

4-3.  Domain  Discretization 

Numerical  solutions  for  the  IB  VP  (4.1)  require  discretization  of  the  spatial 
and  temporal  computational  domain.  Figure  4.1  introduces  a  discretized  x-t 
mesh  convention  where  the  spatial  domain,  of  length  L,  is  partitioned  into  / 
computational  cells  of  uniform  width 

Ax  =  y .  (4-2) 

Each  computational  cell  is  bounded  by  faces  i  —  1/2  and  *  +  1/2  located  at 

Xi- 1/2  =  (*  - 
%i+ 1/2  =  iAx- 

The  center  of  the  ith  computational  cell,  xt ,  is  determined  as 

=  (*  -  ^)Ax-  (4-4) 


(4.3) 
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Figure  4.1  shows  non-uniform  time  steps  which  are  determined,  adaptively,  as 


A  t  = 


c 'cjA % 

qn 

u  max 


(4.5) 


t 


1  .  i-  1  i  i+ 1  .  / 

Figure  4.1.  Discretized  Domain 

Ccji  is  the  Courant  coefficient  satisfying  the  condition 

o  <  Ccfl  <  1 .  (4.6) 

As  Ccfl  approaches  one,  the  time  stepping  scheme  exhibits  maximum  efficiency. 
S]’nax  is  the  maximum  wave  speed  present  throughout  the  spatial  domain  during 
time  level  n.  For  the  time-dependent,  Euler  equations,  a  reliable  approximation 
for  S'xiax .  which  extends  to  multi-dimensional  problems,  (Toro,  1999:221)  is 

Smax  =  max{\u?\  +  a”}.  (4.7) 
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The  conditions  set  forth  in  (4.5)  -  (4.7)  predict  no  wave  present  in  the  solution  of 
all  Riemann  problems  travels  a  distance  greater  than  Ax  during  At .  This 
prediction  correlates  directly  to  solution  stability  during  a  time  step,  and  will  be 
further  discussed  later. 

With  the  x-t  domain  fully  discretized,  solutions  for  (4.1)  ensue.  Results 
are  the  piecewise  constant  distribution  of  cell  averaged  values  assigned  to 
computational  cells  as  (Toro,  1999:215) 

V(x,tn)  =  \J?.  (4.8) 

4- 4-  Discretization  of  the  Euler  Conservation  Equations 

In  order  to  discretize  (4.1)  as  applied  to  the  Euler  conservation  equations, 
discontinuous  solutions  must  be  considered.  In  these  cases,  the  smoothness 
assumption  that  leads  to  the  differential  form  of  the  conservation  equations  no 
longer  holds  true  (Toro,  1999:19).  Therefore,  development  of  a  finite  volume 
discretized  form  of  (4.1)  begins  with  consideration  of  the  integral  form  of  the 
equations  applied  to  a  control  volume  , x2 }  x  [f ,  t2\  in  the  domain  of  interest. 

Integrating  the  system  of  conservation  equations  in  x-t  space  and  using  Green’s 
theorem  yields  the  integral  form  of  the  conservation  equations 

^[udx-F(U)di]=0,  (4.9) 

where  the  line  integration  is  performed,  counter-clockwise,  along  the  boundary  of 
the  domain  (Toro,  1999:62).  This  integral  may  be  expanded  as 
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x2  h  h 

J  U(x,  t2)dx  =  J  l J(x,t±)dx  +  J  F(U(x!,t))dt  -  J  F(U (x2,t))dt  ,  (4.10) 

Xl  X1  \  ti 

which  yields  the  conservative,  discretized  form  of  the  Euler  conservation 

equations  (Toro,  1999:217) 

U“+1  =U"+^[F,--l/2-F,+1/2].  (4.11) 

Simply  stated,  (4.11)  indicates  the  solution  is  determined  by  U"  and  the  net 
flux  through  a  cell  where  Fi_1  represents  the  conserved  variable  flux  entering 
the  cell  through  face  i  —  1/2  and  ~Fl+1//2  is  the  flux  exiting  the  cell  through  face 
7  +  1/2  .  These  numerical  fluxes  are  determined  via  equation  (2.6).  The 
solutions  for  primitive  variables  at  each  cell  face,  required  to  solve  (2.6),  are 
evaluated  by  solving  Riemann  problems  at  the  faces.  Godunov’s  scheme  provides 
a  method  to  determine  all  components  required  for  solutions  of  (4.11)  for  a 
computational  domain. 

4-5.  The  Godunov  Scheme 

Godunov’s  first-order  upwind  method  is  a  conservative  method  of  the  form 
(4.11)  where  intercell  numerical  fluxes  are  computed  using  solutions  of  local 
Riemann  problems.  Spatial  discretization  for  a  given  time  step  is  depicted  in 
Figure  4.2. 
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Figure  4.2.  Domain  Discretization  for  a  Single  Time  Step 


x 


Note  that  each  pair  of  cells  provides  information  to  solve  the  RP,  in  primitive 


variable  form,  at  the  shared  cell  face  i  +  1/2  during  time  step  At  such  that 


W,  +  F(W)X  =  0, 


W(x,  tn)  =  WP 


x  <  0 

x  >  0. 


(4.12) 


Therefore,  exact  solutions  for  (4.12)  at  each  cell  face  provide  the  requisite 
information  to  determine  conserved  variable  fluxes,  via  equation  (2.6).  These 
fluxes,  coupled  with  cell  conserved  variable  values  at  the  beginning  of  a  time 
step,  U”  ,  and  a  specified  time  step  determined  by  (4.5)  facilitate  solutions  for 
cell  conserved  variables  at  the  end  of  a  time  step.  Finally,  these  updated 
cell  conserved  variables  serve  as  initial  conditions  for  the  RPs  during  the 
subsequent  time  step. 

Several  aspects  of  Godunov’s  scheme  require  special  consideration. 
Specifically,  in  order  to  solve  the  RPs  at  x  =  0  and  x  =  L ,  and  update  solutions 
in  cells  1  and  M,  computational  cells  must  exist  to  the  left  of  x  =  0  and  to  the 
right  of  x  =  L  ,  respectively.  Additionally,  determination  of  a  time  step  At  which 
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maximizes  computational  efficiency  and  stability  is  significant.  Finally, 
calculating  primitive  variables  and,  hence,  numerical  fluxes  through  cell  faces 
requires  a  modification  to  the  ERS  sampling  procedure. 

4-6.  Special  Considerations 

4-6.1.  Phantom  Cells 

The  imposition  of  phantom  cells  ensures  that  solutions  at  the 
computational  domain  boundaries  are  consistent  with  conditions  defined  by  the 
physical  model.  When  modeling  shocks,  two  types  of  boundaries  are  considered: 
reflective,  or  symmetric,  and  transmissive. 

If  the  boundary  located  at  x  =  L  represents,  physically,  a  fixed,  reflective, 
impermeable  surface,  then  the  physical  border  is  correctly  modeled  by  a  fictitious 
or  phantom  state  Wp+1  that  is  defined  from  the  known  computational  state  Wp 
such  that 

n  n  „  n  „  n  n  n  /  A  -i  Q\ 

Pi+i  —  Pi  5  ui+ i  —  ~ui  5  Pi+ i  —  Pi-  (4-13) 


Similarly,  a  physically  symmetric  boundary  location  at  x  =  0  requires 


n  n  n  n  n  n 

Po  —  Pl  5  U0  —  ~U1  5  Po  —  Pi  ■ 


(4.14) 


Transmissive  conditions  attempt  to  numerically  reproduce  boundaries  that 
allow  the  physical  passage  of  waves  without  any  effect  on  them.  For  this  case, 
the  primitive  variables  in  phantom  cells  at  x  =  L  and  x  =  0  are  defined  as 

n  n  „  n  „  n  n  n  / A  -i  r\ 

Pi+i  —  Pi  i  ui+ 1  —  ui  5  Pi+i  —  Pi  >  (4-15) 
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and 


n  n  n  n  n  n  n  /A  i 

Po  —  Pi  ,  —  ui  >  Po  —  Pi-  (4-16) 

Although  transmissive  boundaries  work  well  in  a  single  dimension,  attempts  to 
model  their  effects  in  multi-dimensional  problems  produces  false  reflections,  and 
is  an  area  of  significant  ongoing  research  (Toro,  1999:224). 

4  -6. 2.  Time  Step  Stability 

Recall  that  a  stable  time  step  is  approximated  from  the  conditions  defined 
in  (4.5)  -  (4.7).  Now  consider  the  solutions  to  the  Riemann  problems  at  cell  faces 
xi+ 1/2  and  .xi+3/2as  depicted  in  Figure  4.3.  Observe  the  left  traveling  shock, 
resulting  from  the  solution  of  the  RP  at  travels  rapidly  to  the  left  and 

reaches  the  boundary  located  at  27+1/2  before  the  completion  of  time  step  At . 
This  shock  effectively  interferes  with  the  solution  of  the  RP  at  2,/+1/2  and  results 
in  an  erroneous  flux  calculation  that  propagates  outward  toward  other  cells.  As 
incorrect  computations  are  advanced  forward  in  time,  errors  are  compounded 
introducing  numerical  instabilities.  To  remedy  this  problem,  the  time  step  is 
shortened  to  At'  and  interference  between  adjacent  RPs  is  avoided. 

Time  steps  are  easily  shortened  by  modifying  the  value  of  the  CFL 
coefficient.  Recall  the  maximum  wave  speed  is  approximated  by  cell  values  of 
material  velocity  and  speed  of  sound.  For  shock  tests  where  stationary  states 
initially  exist  (it  =  0)  the  maximum  wave  wave  speed  is  dominated  by  the  speed 
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RP(i  +  1/2)  RP(i  +  3/2  ) 


Figure  4.3.  Wave  Patterns  in  Adjacent  Riemann  Problems 

of  sound,  and  is  often  underestimated.  The  result  is  a  time  step  that  is  too  large 
which,  potentially,  introduces  instabilities  early  in  the  computations.  In  order  to 
maximize  stability  and  computational  efficiency,  Toro  recommends  Ccji  =  0.2  for 
the  first  five  time  steps,  and  Cc^  =  0.9  thereafter. 

4  -6. 3.  Sampling 

Chapter  3  presented  the  complete  exact  solution  to  a  general  Riemann 
problem  for  the  Euler  equations.  Recall  that  the  sampling  procedure  requires  a 
coordinate  transformation  from  the  physical  domain  of  the  shock  tube,  x ,  to  the 
computational  domain,  x' ,  of  the  RP.  The  solution  for  primitive  variables  at  a 
given  point  within  the  shock  tube  is  then  determined  by  comparing  x '/ 1  to  the 
wave  speeds  generated  by  the  RP. 

For  Godunov’s  scheme,  sampling  is  performed  at  the  cell  face,  only,  for  the 
special  value  x '/ 1  =  0  .  Recall  the  ERS  determined  that  x '/ 1  <  u*  indicated  the 
sample  point  was  located  left  of  the  contact  wave.  For  Godunov’s  method,  this 
is  analogous  to  u*  >  0.  Similarly,  u*  <  0 indicates  the  sample  point,  or  intercell 
boundary,  is  located  right  of  the  contact  when  applying  Godunov’s  scheme. 
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Finally,  each  of  the  shock  and  rarefaction  wave  patterns  analyzed  during  ERS 
sampling  must  be  taken  into  account  by  Godunov’s  method. 

Godunov’s  scheme  may  be  summarized  as  follows.  Following  application 
of  boundary  conditions,  via  phantom  cells,  a  stable  time  step  is  approximated. 

At  each  cell  face,  a  Riemann  problem  is  solved,  exactly,  in  order  to  determine 
primitive  variables  at  the  face.  The  RP  solutions  allow  determination  of 
conserved  variable  fluxes  through  each  face  from  equation  (2.6).  Finally, 
conserved  variables  in  each  cell  are  updated  by  (4.11)  and  primitive  variables  are 
determined  by  equation  (2.4).  The  scheme  continues,  forward  in  time,  until  the 
desired  solution  time  is  achieved. 

^.7.  Numerical  Results  of  Godunov’s  Method 

Here  the  performance  of  Godunov’s  method  is  analyzed  against  the 
benchmark  shock  tube  test  from  Chapter  3,  consisting  of  a  left  rarefaction, 
contact,  and  right  shock.  The  results  were  determined  using  a  spatial  domain 
discretized  with  I  =  200  computational  cells.  Boundary  conditions  were  set  as 
transmissive;  Toro’s  method  for  CFL  coefficient  selection  was  applied.  Figure  4.4 
shows  a  density  profile  at  t  =  0.25  sec  following  125  time  steps. 

Solutions  are  smeared  at  the  shock  and  contact  discontinuities,  as  well  as 
at  the  head  and  tail  of  the  rarefaction,  which  is  a  common  feature  of  first-order 
methods.  Note  that  the  shock  front  is  smeared  over  four  cells,  while  the  contact 
surface  is  distributed  over  27  cells.  In  general,  contact  waves  are  more 
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Figure  4.4.  ERS  vs.  Godunov  Density  Profile  at  t  =  0.25  sec 


difficult  to  resolve  than  shocks  due  to  the  behavior  of  the  characteristics.  Recall 
the  characteristics  on  either  side  of  a  contact  run  parallel  to  the  discontinuity, 
while  shock  characteristics  run  into  the  shock.  This  compression  mechanism 
effectively  enhances  numerical  shock  resolution.  Lastly,  the  rarefaction  wave, 
which  is  a  smooth  flow  feature,  exhibits  smearing  near  the  head  and  tail,  where  a 
discontinuity  in  the  derivative  exists  (Toro,  1999:227). 

It  is  expected  that  resolution  would  improve  at  each  smearing  location  as 
the  coarse  mesh  of  200  cells  is  refined.  In  fact,  increasing  the  number  of 
computational  cells  by  a  factor  of  five,  from  200  to  1000,  achieves  improved 
resolution  at  solution  and  derivative  discontinuities,  as  depicted  in  Figure  4.5. 
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Figure  4.5.  ERS  vs.  Godunov  Density  Profile  (1000  cells)  at  t  =  0.25  sec 

Here,  the  shock  front  remains  smeared  over  four  cells,  while  the  contact 
discontinuity  is  resolved  over  64  cells.  Because  the  spatial  mesh  refinement  did 
not  increase  smearing  of  the  contact  surface  by  the  same  factor  of  five,  resolution 
of  the  associated  numerical  solutions  is  improved. 

Godunov’s  method,  used  in  conjunction  with  an  exact  Riemann  solver,  is 
shown  to  accurately  predict  the  locations  and  speeds  of  each  wave.  This  feature 
of  the  scheme  is  important  when  modeling  shock  propagation.  Furthermore,  the 
method  guarantees  monotonic  solutions  (Toro,  1999:226)  near  discontinuities, 
which  will  prove  significant  when  implementing  second-order  methods.  Finally, 
practical  computations  involving  Godunov’s  method  require  solutions  to  billions 
of  Riemann  problems,  making  the  iterative  ERS  solution  process  the  most 
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demanding  task  in  the  numerical  method.  In  the  next  section,  improvements  in 


computational  efficiency,  to  include  solution  approximations,  are  examined. 

4 -8.  Adaptive  Riemann  Solvers 

The  computational  effort  required  by  the  iterative  exact  Riemann  solver 
may  not  always  prove  justifiable.  For  this  reason,  several  non-iterative 
approximate  Riemann  solvers  have  been  developed  in  recent  decades.  These 
improved  solution  algorithms  are  based  on  Godunov’s  scheme,  and  offer  the 
benefit  of  improved  efficiency  without  sacrificing  accuracy.  Riemann  problem 
solutions  may  be  estimated  by  one  of  two  methods:  numerical  flux 
approximations  or  state  approximations.  The  latter  form  the  basis  for  modern 
Godunov-type  methods  and  are  examined  here. 

Recall  the  exact  Riemann  solver  requires  finding  the  root  of  Toro’s 
pressure  equation  in  order  to  ultimately  determine  values  for  primitive  variables 
in  the  star  region.  Toro  presents  three  approximations,  based  on  known  left  and 
right  states,  and  an  assumption  of  wave  patterns  in  the  RP,  that  predict  these 
same  variables  with  equal  accuracy  and  greater  efficiency  than  the  exact 
techniques.  Furthermore,  Toro  includes  a  set  of  conditions  that  apply  known 
state  variables  to  adaptively  select  the  most  favorable  value  of  the  determined 
approximations . 

Toro’s  techniques  are  easily  implemented  because  two  of  the  three 
approximate  Riemann  solvers  are  based  upon  the  exact  Riemann  solver. 
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Additionally,  adaptive  selection  of  the  most  favorable  approximation  applies  the 
same  logic  used  to  select  the  initial  pressure  guess  for  Newton’s  method  in  the 
ERS. 

Toro  refers  to  the  three  approximate  Riemann  solvers  as  the  primitive 
variable  Riemann  solver  (PVRS),  two-rarefaction  Riemann  solver  (TRRS)  and 
two-shock  Riemann  solver  (TSRS).  The  logic  sequence  for  selecting  an 
approximate  solver  (Toro,  1999:306)  is  included  in  Figure  4.6,  followed  by 
governing  equations  for  each  solver. 


Figure  4.6.  Logic  Sequence  for  Adaptive  Noniterative  Riemann  Solver 

The  PVRS  relies  on  the  assumption  that  smooth  flows  exist  in  the  region 
of  interest.  This  assumption  yields  a  linear,  hyperbolic  system  of  equations  which 
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is  solved,  exactly,  to  provide  relationships  for  primitive  variables  in  the  star 
region  as 


WPras(WL,WR)  = 


1  1 

P *  =  ~(Pl  +  Pr)  +  ~[(ul  _  ur)(Pl  +  Pr ) (aL  +  %)]? 

Z  o 


u*  =  |(«L  +%)  +  2 


(pl  -  Pr) 


P*  l  —  Pl  + 
P*r  =  Pr  + 


(Pl  +  Pr)(°l  +  °r) 
(^l  ~  u*)(Pl  +  Pr) 

(aL  +  %,) 

(u*  ~  %)(Pl  +  Pr) 

(°l  +  °r) 


(4.17) 


Recall  the  exact  Riemann  solver  pressure  function 

/i(P*,WL)  +  /fl(p„  WR)  +  nR-uL  =0, 


(4.18) 


where  the  left  and  right  functions  are  based  on  the  classification  of  the  outer 
waves  in  the  RP  as 


4(p*,  wc)  = 


(P*  ~  P() 


1 

\~ 

2 


2at 


7-1 


t  Y 
p* 


pe(p*(7  +  l)  +  P^(7-l)) 

7-1 


if  p*  >  p^  (shock), 


Pi 


27 


-1 


(4.19) 


if  p*  <p ( rarefaction ). 


The  TRRS  assumes  both  outer  waves  are  rarefactions  and  substitutes  the 
appropriate  relationships  from  (4.19)  into  (4.18).  The  resulting  closed-form 
solutions  for  primitive  variables  in  the  star  region  (Toro,  1999:301)  are  calculated 
as 
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HW(wL,wR)  = 


P* 


ry  -  \ 

aL  +  aR  77  ( UL  ~  UR ) 

27 


aL  aR 


27 

7-1 


(^)  (^) 

27  27 

Pl  7  Pr  7 
1  1 

«*  =  “(«L  +%)  +  -[/rW-/lW], 


(  \ 


P*  l  —  Pl 


p* 


Pl 


7  V 


P*R  ~  PR 


p* 


PR 


(4.20) 


Similarly,  the  TSRS  classifies  both  outer  waves  as  shocks.  Unfortunately, 
this  approximation  does  not  lead  to  a  closed  form  solution.  Resulting  quadratic 
equations  lead  to  non- uniqueness  of  solutions  and,  in  the  case  of  complex  roots, 
non-existence  of  solutions.  An  acceptable  alternative  is  to  first  estimate  pressure 
as  ppvrs  and  then  apply  the  two-shock  assumption  to  solve  the  pressure  function. 
The  resulting  equations  are 


Wtsrs  (Wl  >  WR ,  p0 ) 


r)  Pl  (Po  )Pl  +  Pr  (Po  )Pr  -  O'r  ~  7 ) 

Pl(Po)  +  Pr(Po) 

u*  =~(UL+  UR )  +  -[(p*  -  Pr)Pr(Po)  -  (P*  -  Pl)Pl(Po)L 


P*L 


P*R 


Pl 


Pr 


P*  [  (7-1) 
Pl  (7  + 1) 
(7  ~  1)  P*  }  x 
(7  + 1)  Pl 
P*  j  (7-1) 
Pr  (7  + 1) 
(7-1 )  P*  ]  t 
.  (7  +  1)  Pr 


(4.21) 
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where 


Po  =  ma,x(0,  pPVRS)  (4.22) 

and 

l 

9t(Po)  = 

The  adaptive  techniques  were  implemented  in  the  computational  model 
and  applied  to  the  benchmark  shock  tube  test.  Parameters  discussed  in  section 
4.7  were  maintained.  A  comparison  of  the  test  results  is  shown  in  Figure  4.7. 


pe[(7  +  i)p0  +(7-1)^ 


(4.23) 


p(kg/m3) 


Figure  4.7.  Comparison  of  Godunov’s  Method  Using  ERS  and  ARS 

While  the  results  generated  by  the  exact  and  approximate  techniques 
appear  graphically  identical,  the  computational  costs  differ  significantly. 
Specifically,  calculations  executed  with  the  approximate  solver,  for  all  five  shock 
tests  presented  by  Toro,  ran,  on  average,  nearly  36%  faster  than  exact 
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computations.  Based  on  this  observation,  I  conclude  that  the  approximate  solver 


is  best  suited  for  all  subsequent  calculations.  Finally,  an  accurate  and  efficient 
first-order  method  is  established  and  second-order  accuracy  is  examined  next. 

4-9.  F.F.  Toro’s  Weighted  Average  Flux  (WAF)  Method 

Here  a  second-order  extension  to  Godunov’s  method  is  explored.  Flux 
calculations  in  Godunov’s  first-order  scheme  require  solutions  to  Riemann 
problems  along  cell  faces  only.  The  weighted  average  flux  method  analyzes  the 
full  structure  of  the  RP  in  order  to  determine  the  total  flux  through  a  cell  face. 
Toro  presents  two  versions  of  his  WAF  method.  The  first  scheme  determines  an 
integral  average  of  the  flux  across  the  full  solution  of  a  local  RP.  The  second 
version  computes  a  weighted  average  state  of  primitive  variables  across  the  local 
RP  and  applies  these  states  to  determine  WAF  (Toro,  1999:492).  The  latter 
version  requires  fewer  flux  calculations  than  the  former,  and  is  implemented  here 
due  to  its  computational  efficiency. 

Weighted  average  flux  is  defined  as  (Toro,  1999:429) 


As  applied  to  the  weighted  average  state  (WAS)  version  of  flux,  (4.24)  may  be 
rewritten 


F*+l/2 


1 

Ax 


Ax/2 

_  /\f 

F(  Wj+i  /  2  (x,  — ) )  dx 
2 

-Ax/2 


(4.25) 
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where  the  WAS  of  primitive  variables  across  a  local  RP  is  expressed 


Wj+l/2 


1 

Ax 


At N 


Ax/ 2 

/  Wm/2(z,^))cfe. 

-Ar/2 


(4.26) 


Figure  4.8  depicts  the  general  wave  structure  of  a  local  Riemann  problem 
with  N  waves,  where  the  flux  weights,  (3k ,  for  k  =  1,  ...,  N+l  are  normalized 
lengths  of  segments  Ak_1Ak  such  that  (3k  =  |^4yfc_1^4fc  |  /  Ax . 


Figure  4.8.  General  Wave  Structure  of  a  Local  Riemann  Problem 

These  weights  can  be  expressed,  in  terms  of  the  wave  speeds,  as 


0k  2 (°k 

AtSk 

ck  ~  ^  5  c0  —  CN+ 1 


1. 


(4.27) 


48 


Here  ck  is  the  local  courant  number  for  wave  k  moving  with  speed  Sk  and 
represents  the  fraction  of  the  cell  width  that  the  kth  wave  traverses  during  a  time 
step  At .  In  expanded  form,  the  weights  are  expressed  as 

A=^(H-ci),  =^(c2  “ci)>  Ab  =^(c3 -G>)>  /?4  =  ^ (!  —  c3 ) •  (4-28) 


For  the  wave  structure  of  Figure  4.8,  the  WAS  integral  (4.26)  becomes 

N+ 1 

E 

k= 1 


Wi+l/2  =  E  A  Wt+l/2 


(4.29) 


where  the  corresponding  states  are  defined  as 

w(1)  =  WL  ,  W(2)  =  W*L  ,  W(3)  =  W*^  ,  W(4)  =  WR  ,  (4.30) 


and 


N+l 

EA=i.  (4.31) 

k= 1 


The  weighted  average  state  of  each  primitive  variable  is  now  fully  defined. 
Weighted  average  flux  is  determined  merely  by  applying  equation  (2.6)  to  the 
WAS  variables. 

To  demonstrate  the  second-order  extension  of  Godunov’s  method,  the  Sod 
shock  test  was  revisited.  Parameters  discussed  in  section  4.7  were  maintained.  A 
comparison  of  test  results  is  shown  in  Figure  4.9.  Inspection  of  the  figure 
indicates  improved  resolution  at  discontinuities  over  Gounov’s  first-order  method. 
Specifically,  smearing  at  the  contact  discontinuity  experienced  a  reduction  from 
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27  cells  to  22  cells.  However,  the  WAF  method  introduced  spurious  oscillations 


at  each  discontinuity,  sometimes  called  Gibb’s  phenomena,  (Toro,  1999:  445) 
which  are  common  side  effects  of  second  order-methods.  Analysis  of  the  WAF 
method  and  Figure  4.8  provides  an  explanation  of  the  fluctuations. 


ptkg/m3) 


Figure  4.9.  WAF  Method  Density  Profile  at  t  =  0.25  sec 

Recall  that  Godunov’s  first-order  method  determines  flux  along  the  cell 
face,  only.  Inspection  of  Figure  4.8  shows,  whether  S2  >  0  or  S2  <  0 ,  that  the 
cell  face  in  the  local  RP  always  lies  in  the  upwind  portion  of  the  star  region. 
With  regard  to  solutions,  upwind  flux  contributions  correlate  to  stability  while 
downwind  contributions  facilitate  increased  accuracy.  Therefore,  because 
Godunov’s  first-order  flux  calculations  contain  only  an  upwind  flux  contribution, 
stable  solutions,  which  may  be  refined  for  accuracy,  are  expected.  The  WAF 
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method  introduces  these  refinements  in  terms  of  downwind  flux  contributions. 
Note  from  Figure  4.8  that  when  S2  >  0  the  WAF  contribution  in  the  star  region 
consists  of  the  upwind  flux  with  weight  J32 ,  and  the  downwind  flux  with 
weight  /?3 .  For  the  case  S2  <  0 ,  the  upwind  flux  weight  in  the  star  region  is  /?3 
while  the  downwind  flux  weight  is  f32 .  In  both  situations,  the  upwind  flux 
contribution  is  larger  than  the  downwind  contribution  and  the  WAF  method  is, 
therefore,  upwind  biased  (Toro,  1999:418).  This  fact  indicates  the  WAF  method 
maintains  the  stability  of  its  first  order  predecessor,  but  benefits  from  increased 
accuracy  due  to  downwind  contributions.  Unfortunately,  these  downwind 
contributions  introduce  oscillations  in  discontinuous  regions  where  steep  gradients 
exist.  E.F.  Toro’s  limited  weighted  average  flux  method  uses  adaptive,  total 
variation  diminishing  techniques  that  effectively  ameliorate  these  oscillations 
while  preserving  stability  and  accuracy. 

4- 10.  E.F.  Toro’s  Limited  Weighted  Average  Flux  (LWAF)  Method 

Total  variation  diminishing  (TVD)  schemes  are  intimately  linked  to 
traditional  artificial  viscosity  methods  as  both  techniques  attempt  to  eliminate  or 
control  false  fluctuations  near  high  gradients.  Whereas  artificial  viscosity 
schemes  introduce  oscillation  reduction  mechanisms  explicitly,  moderation  is 
inherent  in  TVD  methods. 
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Given  a  function  u  =u(x,t ) ,  the  total  variation  of  the  function  at  a  fixed 
time  t  =  tn  is  denoted  TV(u(t )) .  If  it"  =  {«"  j  is  a  mesh  function,  then  the  total 
variation  of  un  is  defined  as  (Toro,  1999:448) 

TV(un)  =  £  - 11  i  •  (4-32) 

i= 1 

A  scheme  is  said  to  be  TVD  if 

TV(un+x)<TV(un),  Vn .  (4.33) 

Finally,  the  above  definitions  produce  the  consequence 

TV(un)  <  TV(un~X )  < ...  <  TV(u°)  .  (4.34) 

Two  important  characteristics  of  TVD  methods  are  convergence  of 
solutions  and  preservation  of  monotonicity.  The  benefit  of  the  former  is  obvious. 
As  for  the  latter,  schemes  that  preserve  monotonic  behavior  of  solutions  predict 
when  data  jw”j  is  monotonic,  the  solution  set  |m”+1|  is  monotone  in  the  same 

manner.  Explicitly,  if  ju”j  is  monotonic  increasing  so  is  |wf+1|  >  an<4  if  \u'i  \  is 
monotonic  decreasing  so  is  j w”+1 1 .  Because  Godunov’s  first-order  method 

guarantees  monotonic  solutions  near  discontinuities,  methods  that  are  TVD  must 
preserve  this  behavior  in  the  presence  of  strong  flows. 

In  order  to  limit  downwind  flux  contributions  and  ensure  monotonic 
behavior  near  discontinuities,  Toro  presents  four  weight  limiter  functions  that 
effectively  regulate  flux  weights  within  downwind  regions  of  local  Riemann 
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problems.  These  weight  limiter  functions,  ^(r,  c  )  ,  are  derived  directly  from 
corresponding  flux  limiter  functions  y/(r)  and  are  related  as  (Toro,  1999:499) 

<f)(r, \c\ )  =  1  -  (1  - \c\)y/(r)sign{c) .  (4.35) 


Here,  c  is  the  local  courant  number  given  by  equation  (4.27)  and  r  represents  a 
flow  parameter  that  compares  upwind  and  local  density  gradients  as 


(4.36) 


These  equations  represent  scalar  relationships  applied  to  a  single  wave.  In 
practice,  solutions  of  the  Euler  conservation  equations  in  a  single  dimension 
produce  three  waves  resulting  from  the  three  conservation  laws.  Therefore, 

(4.35)  and  (4.36)  are,  in  fact,  applied  across  each  wave  in  the  solution  of  the 
Riemann  problem. 

Before  introducing  the  governing  equations  associated  with  flux  limiters, 
an  illustrative  examination  of  flow  parameter  calculation  is  presented  in  Figure 
4.10.  For  each  case,  the  local  density  gradient  at  the  i  +1/2  boundary  across  the 
kth  wave  is  calculated 

Api+i/2,h  =  Pi+l/2,k+l  —  Pi+l/2,k  •  (4.37) 


For  the  case  where  ci+1/2 ,k  >  0,  the  upwind  location  lies  left  of  the  local  position 

and  the  upwind  density  gradient  is  determined  at  the  i  -  1/2  boundary  across  the 
kth  wave  as 
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APi-i/2,k~  Pi— 1/2,  k+1  Pi—1  /2,k  ' 


(4.38) 


Conversely,  when  ci+1/2,k  <  0»  the  upwind  site  lies  right  of  the  local  location  and 

the  upwind  density  gradient  is  calculated  at  the  i  +  3/2  boundary  across  the  kth 
wave  as 

APi+3/2,k  =  Pi+3/2,k+l  ~  Pi+3/2,k  •  (4.39) 


Determination  of  the  flow  parameter  introduces  a  side  effect  of  the  LWAF 
method  that  requires  special  consideration.  For  the  case  where  ci+1/2)  k  >  0  and 

the  local  boundary  i  +  1/2  lies  at  the  left  edge  of  the  physical  domain  at  x  =  0, 

ci  +  n.t>0  ci  +  i/2k<0 

Wave  i  —  1/2,  k  Wave  i  + 1/2,  k  Wave  i  + 1/2,  k  Wave  i  +  3/2,  k 


i  — 1/2  i  + 1/2  i  +  3/2  i  — 1/2  i  + 1/2  I  +  3/2 

Figure  4.10.  Determination  of  Flow  Parameter 

the  upwind  density  gradient  across  boundary  i  —  1/2  requires  solution  of  a 
Riemann  problem  via  phantom  cells  0  and  -1.  Although  phantom  cell  0  exists, 
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cell  -1  must  be  added.  Similarly,  when  ci+]y2)fc  <  Oand  the  local  boundary  at  i  + 

1/2  lies  at  the  right  edge  of  the  physical  domain  at  x  =  L,  the  upwind  density 
gradient  across  boundary  z  +  3/2  requires  solution  of  a  Riemann  problem  using 
phantom  cells  1+1  and  1+2.  Here,  phantom  cell  1+2  is  added. 

To  correctly  model  physically  reflective  or  symmetric  boundary  conditions 
at  x  =  0  and  x  =  L  these  added  phantom  cells  are  defined  from  known 
computational  states  as  (Toro,  1999:491) 


n  n  n  n  n  n 

P- 1  =  P2  ,  U-1  =  -«2  5  P-1  =  P2 

n  n  „  n  nn  n  n 

Pl+ 2  —  Pi- 1  >  UI+ 2  —  ~UI- 1  >  Pl+ 2  —  P/-1 


(4.41) 


Likewise,  transmissive  boundaries  are  given  by 


n  n  „  n  _  „  n  _  „n 

P-1  —  P2  i  U-1  —  >  P-1  —  P2 


„n  „  n  „  n 

Pi +2  —  Pi- i  >  u/+2  —  M/-i 


P/+2  =  Pi- i 


(4.42) 


Toro  presents  four  weight  limiters  identified  as  Super  A,  Van  Leer  A,  Van 
Albada  A,  and  Mini  A  and  associates  flux  limiters  Super  Bee,  Van  Leer,  Van 
Albada,  and  Mini  Bee  (Toro,  1999:469)  with  the  respective  weight  limiters. 

These  flux  limiters  are  defined  as 


Ab  0) 


0,  r  <  0, 

2 r,  0  <  r  <  1  /  2, 

1,  1  /  2  <  r  <  1, 

r,  1  <  r  <  2, 

2,  r  >  2, 


(4.43) 


55 


0,  r  <  0, 

^vl(r)  =  -  2  r  (4.44) 

7~  ’  r  >  °> 

[1  +  r 

0,  r  <  0, 

VVa(r)  =  '  r(l  +  r)  (4.45) 

1  +  r2 

0,  r  <  0, 

^m6(r)  =  -r>  0  <  r  <  1,  (4.46) 

1,  r  >  1. 

and  are  plotted  in  Figure  4.11. 


0.5  1  1.5  2 

Figure  4.11.  Flux  Limiter  Comparison 


The  flow  parameter,  coupled  with  flux  and  weight  limiter  functions, 
provides  adaptive  adjustments  to  flux  weight  calculations  based  upon  local 
conditions  resulting  from  Riemann  problem  solutions.  The  flow  parameter 
effectively  measures  the  smoothness  of  solutions.  When  the  upwind  and  local 
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density  gradients  are  comparable,  r  approaches  one,  flows  are  characterized  as 


smooth,  and  downwind  flux  contributions  are  included  to  guarantee  second-order 
accuracy  of  solutions  (Toro,  1999:457).  However,  as  r  >>  1,  flows  become 
discontinuous,  imposing  adaptive  adjustment  to  flux  calculations.  Inspection  of 
flux  limiter  functions  (4.43)  -  (4.46)  and  Figure  4.11  shows  fluxes  are  regulated  as 
r  becomes  large.  Specifically,  when  r  exceeds  one,  upwind  flux  contributions  are 
maximized  to  ensure  stability,  while  downwind  flux  contributions  are  limited  in 
order  to  restrict  the  downwind  flux  components  that  produce  false  fluctuations 
near  discontinuous  regions. 

With  TVD  methods  fully  introduced,  the  weighted  average  state  of  each 
primitive  variable  must  be  determined.  A  TVD  scheme,  analogous  to  equation 
(4.29)  is  given  as 

_  N+l 

W.+1/2  =  ^  Pk,LWAF  W<‘>1/2  (4.47) 

k=  1 


where 

L  WAF  =  -  +  ^1^5'n(cl))’ 

& ,lwaf  =  \{(kSign(c 2)  -  (faSignfa)) , 

1  (4.48) 

Pz,LWAF  =  -(03%n(c3)  -  (hSignfo)), 

Pa,LWAF  =  2  (1  ■  03^fi,n(c3))- 
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Finally,  limited  weighted  average  flux  is  determined  merely  by  applying  equation 
(2.6)  to  the  LWAS  variables. 

Application  of  the  flux  and  weight  limiters  to  the  benchmark  shock  tube 
test  produced  results  depicted  in  Figure  4.12,  Figure  4.13,  and  Figure  4.14. 

Figure  4.12  indicates  that  each  flux  limiter  effectively  eliminates  spurious 
oscillations  near  the  shock  and  contact  surface  discontinuities.  Figure  4.13  and 
Figure  4.14  reveal  detailed  behavior  of  each  limiter  function  in  the  two  regions  of 
interest.  Figure  4.13  indicates  that  the  Super  Bee  flux  limiter  function  performs 
best  at  the  contact  surface,  while  the  Van  Leer  function  provides  more 
graphically  consistent  results,  compared  to  exact  solutions,  than  the  remaining 
three  functions  at  this  discontinuity.  Inspection  of  Figure  4.14  demonstrates 
sporadic  behavior  of  the  Super  Bee  function  at  the  shock  front,  while  the  Van 
Leer  limiter  best  eliminates  oscillations  at  the  discontinuity.  Therefore,  the  Van 
Leer  flux  limiter  is  presumed  best  for  shock  modeling  and  will  be  implemented  for 
all  subsequent  calculations. 
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Figure  4.12.  Flux  Limiter  Sod  Test  Density  Comparison 
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Figure  4.13.  Contact  Surface  Flux  Limiter  Density  Comparison 
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Figure  4.14.  Shock  Front  Flux  Limiter  Density  Comparison 


4-11-  Verification  of  Second- Order  Accuracy 


The  primary  objective  of  numerical  techniques  is  to  determine  accurate 
approximations  to  exact  solutions.  The  order  of  accuracy  of  an  implemented 
scheme  measures  the  amount  by  which  exact  solutions  differ  from  discrete 
numerical  solutions  such  that  (Slater,  2000) 

E  =  \f(h)-fe„cl\  =  Ch“.  (4.49) 


E  represents  absolute  error;  f(h)  and  fexact  symbolize  approximate  and  exact 
solutions,  respectively;  C  is  a  constant;  and  h  is  a  measure  of  grid  spacing  with 


order  of  accuracy  p  where 


hp 


(b  —  a 


\p 


n 


(4.50) 
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Here,  a  and  b  represent  lower  and  upper  domain  limits,  while  n  denotes  the 


number  of  subintervals  used  to  discretize  [a,  b]. 

E  _  |/(^)  —  fexact  | 

hP  ~  hP 


It  follows  from  (4.49)  that 

=  C, 


(4.51) 


which  predicts  that  finite  volume  error  calculations  for  p  =2  using  varying  mesh 
sizes  are  expected  to  yield  a  constant,  for  a  second-order  accurate  scheme. 

Confirmation  of  second-order  accurate  solutions  was  examined  by 
comparing  exact  solutions  from  the  benchmark  shock  test  against  solutions 
calculated  using  the  TVD  LWAF  Godunov  scheme  with  several  mesh  sizes. 
Sample  points  were  selected  from  each  of  the  four  Riemann  problem  states  and 
analyzed  for  error  and  order  accuracy.  Table  4.1  depicts  this  comparison  in  the 
left  star  state. 


n 

E 

E/h 

E/h2 

25 

0.00133348 

0.033337 

0.833425 

50 

0.00232385 

0.116192 

5.80962 

200 

0.000166184 

0.0332369 

6.64738 

400 

0.000031504 

0.0126019 

5.04076 

800 

4.78943  x  10"6 

0.00383155 

3.06524 

1600 

1.33387  x  10‘6 

0.00213420 

3.41472 

3200 

5.98817  x  10‘7 

0.00191621 

6.13188 

6400 

4.07474  x  10‘7 

0.00260783 

16.6901 

12800 

2.64006  x  10'7 

0.00337928 

43.2547 

Table  4.1.  Error  Comparison  for  Density  in  Left  Star  State 


If  the  scheme  is,  in  fact,  second-order  accurate,  the  results  in  column  four  of 
Table  4.1  are  expected  to  converge  to  a  constant.  Similarly,  a  first-order 
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accurate  scheme  would  produce  convergent  results  in  column  three  of  the  table. 


Column  two  of  the  table  indicates  that  solutions  begin  to  converge  as  the  mesh 
size  is  refined  from  500  to  1000  cells.  However,  at  first  glance,  neither  column 
three  nor  column  four  appears  to  converge  to  a  constant  along  this  interval,  and 
the  order  of  accuracy  of  the  scheme  cannot  be  definitively  stated. 

Because  the  magnitudes  of  the  values  in  the  respective  columns  differ  by 
several  orders,  the  error  was  further  analyzed  using  a  log-log  plot.  Specifically, 
taking  the  log  of  both  sides  of  equation  (4.49)  yields 

log  E  =  log  (7  +  plog  h  (4-52) 

which  indicates  that  the  logarithmic  error  varies  linearly  with  mesh  size,  h,  and 
slope  p.  Figure  4.15  shows  data  sets  generated  from  (4.52),  using  p  =  1,  p  =  2, 
and  mesh  sizes  from  Table  4.1,  as  compared  to  error  calculated  from  the  left  star 
state  of  the  computational  model  with  the  same  parameters.  The  figure  indicates 
error  data  generated  using  the  developed  one-dimensional  code  more  closely 
resemble  data  produced  from  (4.52)  with  p  =  2.  The  same  analysis  was 
conducted  using  data  from  the  right  star  state,  and  produced  similar  results. 

Because  the  generated  results  show  better  graphical  correlation  to  second- 
order  accuracy,  as  compared  to  first-order  accuracy,  the  TVD  LWAF  Godunov 
scheme,  empirically,  exhibits  second-order  accuracy. 
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Figure  4.15.  Log-Log  Plot  of  Error  vs.  Number  of  Cells  in  Left  Star  State 
4’ 12.  Validation  of  the  One- Dimensional  Shock  Code 
4-12.1.  Sod  Test  in  One  Dimension 

During  each  step  of  code  development,  calculated  data  were  compared 
against  the  Sod  shock  test  data  computed  with  an  exact  Riemann  solver. 
Graphical  evaluation  of  these  data  sets  verified  proper  implementation  of  each 
introduced  numerical  technique.  In  addition  to  the  Sod  test,  generated  data  from 
four  diverse  shock  tests,  discussed  in  Chapter  3,  were  compared  graphically  to 
results  produced  by  the  ERS.  Numerical  results  proved  consistent  with  published 
data  for  all  tests.  Finally,  in  order  to  guarantee  accuracy,  validation  of  the 
developed  code  against  experimental  data  is  required. 
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4-12.2.  Army  Research  Laboratory  (ARL)  57cm  Shock  Test 


In  1996,  the  Army  Research  Laboratory  (ARL)  used  measurements 


obtained  in  their  57cm  shock  tube  facility  to  validate  a  second-order  accurate 
single  dimensional  shock  code.  The  shock  tube  was  100m  long  with  a  diameter  of 
57cm  and  contained  a  driver  region  0.91m  long  (Schraml,  1996;  Wittig,  1999:72). 
The  following  conditions  existed  at  the  beginning  of  the  experiment. 

[WL  =  (4.486 kg/m3,  0.0 m/sec,  379.2  x  103  N/m2)  x  <  0, 


W(x,  0)  = 


WR  =  (1.208 kg/m3,  O.Om/sec,  102.1  x  103  N/m2)  x  >  0. 


(4.53) 


An  overpressure  history  was  collected  from  a  measurement  station  located 
31.48m  from  the  initial  discontinuity.  Figure  4.16  portrays  overpressure  history 
data  determined  experimentally  and  computationally  by  ARL.  The  shock  front 
arrived  at  the  experimental  measurement  station  at  66.0ms  with  a  measured 
overpressure  of  66.3kPa  (Schraml,  1996;  Wittig,  1999:73). 

Using  the  developed  LWAF  one-dimensional  code  to  accurately  model  the 
ARL  experiment  and,  ultimately,  validate  the  overpressure  at  the  shock  front 
requires,  first,  time  of  arrival  calibration.  To  verify  correct  time  of  arrival  of  the 
shock  front,  a  symmetry  boundary  was  placed  at  x  =  0m  while  a  transmissive 
boundary  was  positioned  at  x  =  50m.  This  far  right  boundary  location  reduced 
execution  time  without  affecting  results  at  the  computational  measurement 
station  at  x  =  31.48m  (Wittig,  1999:74).  The  computational  domain  was 
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Figure  4.16.  ARL  Overpressure  History  at  31.48m 

discretized  into  5000  cells  during  the  time  interval  [0,  100]  msec,  a  CFL 
coefficient  of  0.9  was  applied,  and  the  van  Leer  flux  limiter  function  used. 

Code  execution  using  these  initial  conditions  generated  a  later  time  of 
arrival,  by  0.61ms  (0.9%  relative  error),  and  a  lower  static  overpressure,  by 
6.6kPa  (10.0%  relative  error),  than  observed  during  the  ARL  experiment.  To 
resolve  these  discrepancies,  the  code  was  re-executed  using  varying  mesh  sizes 
and  CFL  coefficients.  Table  4.2  summarizes  the  results  of  mesh  refinement,  using 
a  CFL  coefficient  of  0.9,  for  all  computations. 

Because  mesh  refinement  did  not  significantly  alter  shock  time  of  arrival, 
calculations  were  repeated  using  several  CFL  coefficients  and  a  mesh  size  of  5000 
cells,  to  maximize  computational  efficiency.  Table  4.3  summarizes  the  results  of 
these  calculations. 
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n 

Time  steps 
( CFL  =  0.9) 

Time  of 
Arrival  (msec) 

Max  Static 
Overpressure 
(kPa) 

5000 

5623 

66.61 

59.7 

6000 

6750 

66.62 

59.3 

7000 

7875 

66.63 

59.9 

8000 

8999 

66.69 

59.7 

9000 

10123 

66.70 

59.4 

10000 

11250 

66.71 

59.9 

Table  4.2.  Mesh  Refinement  Effects  on  Time  of  Arrival  and  Static  Overpressure 


CFL 

Time  steps 
(CFL  =  0.9) 

Time  of 
Arrival  (msec) 

Max  Static 
Overpressure 
(kPa) 

0.9 

5623 

66.61 

59.7 

0.8 

6310 

66.54 

58.8 

0.6 

8410 

66.48 

58.8 

0.4 

12609 

66.44 

58.7 

0.2 

25208 

66.45 

58.6 

0.1 

50412 

66.49 

58.6 

Table  4.3.  CFL  Refinement  Effects  on  Time  of  Arrival  and  Static  Overpressure 

As  with  mesh  refinement,  CFL  coefficient  modification  did  not  produce 
significant  changes  to  time  of  arrival  or  static  overpressure.  Because  the 
objectives  of  these  refinements  were  to  decrease  time  of  arrival  and  increase 
overpressure,  data  from  the  above  tables  were  examined  to  predict  the  most 
favorable  combination  of  mesh  size  and  CFL  coefficient.  Based  on  this  analysis, 
a  final  computation  was  run  using  a  mesh  of  7000  cells  with  a  CFL  coefficient  of 
0.4.  Again,  solutions  failed  to  improve.  Finally,  in  order  to  minimize  the 
possibility  of  false  reflections  at  the  computationally  transmissive  boundary, 
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several  simulations  were  completed  by  modifying  the  location  of  the  far  right 
boundary  from  60m  to  100m.  Solutions  did  not  improve,  and  I  concluded  that 
the  discrepancies  resulted  from  the  simplifying  assumptions,  inherent  in  the  Euler 
conservation  equations,  which  neglect  the  physical  effects  of  viscosity  and  heat 
transfer. 

In  the  end,  achieving  the  correct  arrival  time  required  increasing  the  initial 

conditions  in  the  driver  region  by  five  percent  to 

WL  =  (4.710kg/m3,  O.Om/sec,  398.1  x  103  N/m2)  x  <  0, 
W(x,0)=  „  ,  9  (4.54) 

WR  =  (1.208 kg/m3,  O.Om/sec,  102.1xl03  N/m2)  x  >  0. 

Figure  4.17  illustrates  overpressure  history  data  at  the  computational 
measurement  station.  As  depicted  in  Figure  4.18,  compared  to  ARL 
computational  results,  the  LWAF  generated  data  proves  smooth  at  the  shock 
front  (a  result  of  application  of  the  TVD  flux  and  weight  limiter  functions). 

The  developed  LWAF  single  dimensional  code  predicted  shock  arrival  at 
the  computational  measurement  station  after  66.03ms,  and  a  peak  overpressure 
at  the  shock  front  of  62.4kPa.  Because  the  predicted  overpressure  agreed  with 
the  experimental  overpressure  within  5.9%,  I  concluded  the  developed 
computational  model  correctly  models  shocks.  Confident  that  the  implemented 
single-dimensional  techniques  are  sound,  the  numerical  procedures  are  next 
extended  into  two  and  three  dimensions. 
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Figure  4.17.  LWAF  Overpressure  History  at  31.48m 


Time  (ms) 


Figure  4.18.  Overpressure  History  Comparison  at  31.48m 


5.  Code  Implementation  in  Multiple  Dimensions 

5.1.  Overview 

The  fully  developed,  single-dimensional,  second-order  accurate  LWAF 
solver  is  applied  to  two  and  three-dimensional  schemes  for  non-linear  hyperbolic 
systems  of  equations.  Finite  volume  techniques  are  introduced  and  applied  to  the 
time-dependent  Euler  conservation  equations  in  multiple  dimensions.  Special 
considerations  and  improvements  to  computational  efficiency  are  addressed. 
Finally,  the  developed  multi-dimensional  computational  models  are  validated 
against  the  one-dimensional  benchmark  shock  tube  test,  as  well  as  mild 
cylindrical  and  spherical  shock  tests,  and  compared  to  published  results. 

5.2.  Implementation  in  Two  Dimensions 
5.2.1.  Initial  Boundary  Value  Problem 

The  two-dimensional  initial  boundary  value  problem  for  non-linear 
systems  of  hyperbolic  conservation  laws,  assuming  uniform  flow  in  the  ^-direction 
is 

PDEs:  Uf+F(l4+  0(^=0, 

m  (5.1) 

ICs:  V(x,y,0)  =  lJ{0>{x,y). 

Here  U^(a;,  y)is  a  piecewise  continuous  distribution  of  initial  data  defined  over 
the  spatial  domain  [0,  Lx]  x  [0,  L  ]  at  t  =  0.  Finally,  in  two  dimensions,  boundary 
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conditions  are  described  along  a  line  at  each  of  the  four  boundaries  using 
symmetric  and  transmissive  boundaries,  as  described  in  the  one-dimensional  case. 

5.2.2.  Discretization  of  the  Euler  Conservation  Equations 

Discretization  of  the  spatial  and  temporal  computational  domain  [0,  LJ  x 
[0,  L  ]  x  [0,  T]  requires  extension  of  the  single-dimensional  mesh  in  the  y- 
direction.  The  spatial  mesh  is  partitioned  into  /  x  J  computational  cells  of 
uniform  area  Ax  x  Ay  such  that 


Ax  =  ^l, 

I 

Ay  =  — . 
y  J 


(5.2) 


Each  computational  cell  is  bounded  by  faces  i  -  1/2  and  i  +  1/2  in  the  re¬ 
direction,  and  faces  j  -  1/2  and  j  +  1/2  in  the  y-direction  where 

**- 1/2  =  (*  - 
1/2  iAx, 

Vj- 1/2  =  (/  -  !)^ 

^•+1/2  =  jav- 


The  center  of  each  computational  cell  (i,  j)  is  determined  as 

=  ((*  -  bAx,  U  -  bAv)  ■ 


(5.4) 


As  in  the  single-dimensional  case,  the  temporal  domain  is  separated  into  n  non- 
uniform  time  steps  which  are  determined,  adaptively,  as 
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At  =  min 


(5.5) 


CcflAx 

' Qn\ 

,i  J  max 


CcfAy ' 

(S")ra«, 


Similarly,  ^maa.  are  approximations  for  maximum  wave  speeds  in  the 

x  and  y  directions,  respectively,  during  time  level  n.  and  are  based  upon  the 
material  and  sound  speeds  where 


('-h  )max  max  jj'fq  ^  |  T  | , 
(^)max  =  max{|<j|  +  ai,j} 


The  conditions  established  in  (5.5)  and  (5.6),  coupled  with  the  CFL  coefficient 
restriction  defined  in  Chapter  4,  predict  waves  in  the  solution  of  Riemann 
problems  that  remain  spatially  confined  within  cell  limits  defined  in  (5.3)  for 
every  time  step.  This  prediction  correlates  directly  to  solution  stability  in 
multiple  dimensions. 

5.2.3.  Dimension  Splitting 


The  dimension  splitting  technique  permits  resolving  the  two  dimensional 


IBVP  (5.1)  as  two  single-dimensional  IBVPs.  A  simple  version  of  this  approach 
replaces  (5.1)  by  the  sequence  of  IVPs  (Toro,  1999:541) 


PDEs  :  U,  +  F(U)X  =  0, 


ICs  : 


un. 


^ljn+l/2 


hJ 


(5.7) 


and 
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(5.8) 


PDEs  :  U,  +  G(U)„  =  0, 

-^>Un+1 

ICs:  V?t‘/2  ' 

hJ 

The  solution  U”+1//2  is  obtained  via  (5.7)  by  solving  J  one-dimensional  problems 
in  the  ^-direction.  This  solution  provides  initial  conditions  to  subsequently  solve 
(5.8)  using  I  one-dimensional  problems  in  the  ^-direction  in  order  to,  ultimately, 
determine  U™^1 .  In  operator  notation  (5.7)  -  (5.8)  may  be  written 

vn+l  =  Y(At)x(At)^jn^  (5.9) 

where  X^At^  and  F^'*  represent  approximate  time  step  solution  operators  to  the 
IB  VP.  There  is  no  specific  purpose  for  applying  the  operators  in  the  order 
described.  Accordingly,  an  equivalent  scheme  is 

U^1  =  X{At)Y{At\ UV)  (5.10) 

It  can  be  shown  that  (5.9)  and  (5.10)  are  first-order  accurate  during  each 
time  step,  given  that  the  individual  operators  X^At)  and  F1  At’  arc  at  least  first- 
order  accurate  (Toro,  1999:542).  Because  the  desired  solution  must  be,  spatially, 
second-order  accurate,  Toro  presents  two  schemes  that  satisfy  the  accuracy 
requirement  every  other  time  step  as 

u”+2  =  x^Y{At)Y{At)X{At\ uy  (5.11) 

and 

U™+2  =  Y^At)X^At)X^At)Y^At\\]lj) .  (5.12) 
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In  practice,  (5.11)  is  implemented  over  the  interval  [tn,tn+2]  based  upon 
the  following  method.  Given  initial  conditions  ,  boundary  conditions  are 
established  as  symmetric  or  transmissive  and  a  stable  time  step  is  approximated 
via  (5.5).  The  first  of  two  time  steps  begins  with  operating  on  TJX . 

Specifically,  in  the  x-direction,  J  one-dimensional  Riemann  problems  are  solved, 
fluxes  are  calculated  at  each  computational  boundary,  and  temporary  variable 
solutions  are  updated  in  each  cell  as  XJij  .  Next,  boundary  conditions  are 
updated  along  the  lines  y  =  0  and  y  =  Ly  and  Yl^r> operates  on  by  solving 
I  one-dimensional  problems  in  the  ^-direction,  as  described  above.  At  the 
conclusion  of  this  first  time  step,  solutions  are  analyzed  against  the  CFL  stability 
criteria.  If  unstable  solutions  exist,  interim  solutions  are  discarded,  the  CFL 
coefficient  and  time  step  are  re-calculated,  and  the  time  step  is  reset.  However, 
provided  stable  solutions  exist,  the  second  time  step  begins. 

Procedures  in  the  second  time  step  are  executed  identically  to  those  in  the 
first  step,  with  one  exception.  Now  operates  on  ,  yielding  the 

temporary  solution  ,  upon  which  X'Al-  subsequently  operates.  This  reversal 
of  operations  effectively  eliminates  the  y-directional  biased  solution  U y"1  present 
at  the  end  of  the  first  time  step  and  produces  the  desired  second-order  accurate 
solution  U"+2at  the  conclusion  of  the  second  time  step.  Subject  to  solution 
stability  verification,  simulation  time  is  advanced  by  2 At .  Finally,  this  second- 
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order  accurate  scheme  is  executed,  forward  in  time,  until  the  desired  solution 
time  is  reached. 

5.2.4 •  Special  Considerations 

The  dimension  splitting  technique  facilitates  solving  the  Euler 

conservation  equations,  in  the  ^-direction,  as 

pu 

pu 2  +  p 

+  =0,  (5.13) 

E(u  +  p ) 

pvu 


and  in  the  ^-direction 


p 

pv 

pv 

to 

+ 

E 

+ 

E{v  +  p) 

pu 

t 

puv 

(5.14) 


Note  the  first,  second,  and  third  equations  in  system  (5.13)  appear  identical  to 
the  pure  one-dimensional  problem  in  the  x-direction.  In  this  system,  u  represents 
the  normal  component  of  material  velocity,  while  v  represents  the  tangential,  or 
shear,  component.  Similarly,  in  system  (5.14),  v  denotes  the  normal  velocity 
component,  while  u  symbolizes  the  tangential  component.  Consideration  of 
tangential  velocities  introduces  corresponding  tangential  momentum  terms  and 
fluxes  which  appear  in  the  fourth  equations  of  both  systems.  Furthermore, 
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inspection  of  the  systems  demonstrates  these  shear  velocities  contribute  to  total 
energy  density  determination  as 

E  =  -p(u2  +  v2  +  w2)  +  ^ — .  (5.15) 

2  7  —  1 

Solutions  to  (5.13)  and  (5.14)  that  apply  TVD  flux  limiter  functions 
produce  the  dissipative  wave  structure  of  Figure  5.1. 


Figure  5.1.  Two-Dimensional  TVD  Flux  Limiter  Wave  Structure 

This  arrangement  dictates  normal  and  tangential  flux  calculations,  which  are 
completed  separately.  Recall  from  Chapter  4  that,  in  a  single  dimension, 
conserved  variables  fluxes  are  calculated  by  determining  the  weighted  average 
state  of  primitive  variables.  In  the  two-dimensional  case,  the  weighted  average 
states  of  density,  normal  velocity,  and  pressure  are  determined  using  equations 
(4.43)  and  (4.44).  These  primitive  variable  weighted  average  states  permit 
computation  of  fluxes  for  mass,  normal  momentum,  and  energy. 
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Calculation  of  the  tangential  velocity  WAS  and  the  corresponding 
tangential  momentum  flux  and  energy  flux  requires  special  treatment.  To 
determine  the  former  requires  application  of  the  equations  (Toro,  1999:550) 

V  T  ,i+l  /  2  =  P5,LWAfVt,1+  l^6,LWAF^T,4  (5.16) 

where  VT  represents  the  tangential  velocity  and  the  TVD  weights  are 

Pb,LWAF  =  (1  +  04*S^(c4))  (5.17) 


and 

Pq,LWAF  =  (1  —  faSignfa))  •  (5.18) 


The  eigenvalue  structure  of  solutions  to  (5.13)  and  (5.14)  conveniently  dictates 
c4  =  C2  (Toro,  1999:548).  However,  determination  of  </>4  is  not  as  trivial.  Unlike 
weight  limiter  functions  0l5  02,  and  03  ,  which  depend  on  the  density  gradient 

across  the  local  and  upwind  waves,  04(74)  is  calculated  from  tangential  velocity 
gradients  AVT  as 


ri+l/2,4 


/^T[i+l/2-(Sign(c4)),A] 

AVT[i+ 1/2,4] 


(5.19) 


Calculation  of  all  components  required  to  solve  systems  (5.13)  and  (5.14)  is 
complete. 
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5.2.5.  Sod  Test  in  Two  Dimensions 


Because  the  benchmark  shock  tube  test  consists  of  flows  in  a  single 
direction,  correct  modeling  of  planar  flows  in  two  dimensions  is  readily  verified. 
Accordingly,  the  two  dimensional  Sod  test  was  executed  over  the  square  domain 
[0,  lm]  x  [0,  lm]  consisting  of  100  cells  in  the  x  and  y  directions  and  a 
diaphragm,  positioned  along  the  line  x  =  0m.  Initial  conditions  were  established 


as 


W(x,y,0)  = 


WL  =  (1.000  kg/m3, 0.0 m/sec,  0.0m/sec,1.0N/m2)  x  <  0.5m, 
W R  =  (0.125 kg/m3, 0.Om/sec,  O.Om/sec,  0.1N/m2)  x  >  0.5m. 


(5.20) 


and  solutions  developed  over  the  time  interval  t  =  [0,  0.25  sec].  The  density 
profile  at  t  =  0.25  sec  is  depicted  below. 


P 


Figure  5.2.  Two-Dimensional  Sod  Density  Profile  at  t  =  0.25  sec 
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To  ensure  planar  flows  modeled  in  two  dimensions  are  consistent  with  the 


one-dimensional  Sod  test  density  data  from  the  two-dimensional  test  was 
extracted  from  the  first  set  of  cells  in  the  ^/-direction.  These  results  are  compared 
to  the  exact  solution,  as  well  as  the  second-order  accurate  solution  obtained  by 
the  LWAF  method  in  a  single  dimension,  in  Figure  5.3.  The  second-order 
accurate  solutions  appear  graphically  identical  for  both  dimensions. 


p  (kg/m3) 


Figure  5.3.  One  and  Two-Dimensional  Sod  Density  Profiles  at  t  =  0.25  sec 
5.2.6.  Mild  Cylindrical  Shock  Test 

To  verify  that  the  developed  code  correctly  models  flows  in  two 
dimensions,  the  Euler  equations  were  solved  over  the  square  domain  [0,  2m]  x  [0, 
2m]  discretized  as  a  computational  mesh  of  100  cells  in  the  x  and  y  directions. 
The  initial  conditions  consisted  of  a  circular  region  with  radius  0.4m,  centered  at 
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(1,1),  with  primitive  variables  of  and  WR  from  (5.20)  representing  states 
inside  and  outside  the  circle,  respectively  (Toro,  1999:583).  Figure  5.4  shows  the 
density  profile  at  t  =  0.25  sec.  Inspection  of  the  profile  indicates,  as  expected,  a 
solution  which  exhibits  a  circular  shock,  contact  surface,  and  rarefaction. 


Figure  5.4.  Cylindrical  Explosion  Density  Profile  at  t  =  0.25  sec 

The  density  profile  indicates  the  presence  of  potential  numerical  artifacts 
in  the  region  between  the  rarefaction  and  contact  discontinuity.  Accordingly, 
graphical  resolution,  in  addition  to  circular  symmetry,  of  solutions  may  be 
further  examined  through  the  contour  plot  in  Figure  5.5.  Again,  the  artifacts 
appear  in  the  region  of  concern,  where  the  solutions  appear  to  lack  resolution. 
Doubling  the  cells  in  both  directions  from  100  to  200  reduced  the  oscillations  in 
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the  region  of  concern,  as  shown  in  Figure  5.6.  As  expected,  the  graphical  output 
shows  fewer  oscillations  and  better  circular  symmetry  with  mesh  refinement. 


Figure  5.5.  Cylindrical  Explosion  Density  Contours  at  t  =  0.25  sec 

The  same  test  was  conducted  using  a  mesh  of  400  cells  in  each  direction. 
Doubling  the  cells  again  did  not  appear  to  significantly  change  solutions,  either 
graphically  or  numerically.  Additionally,  increasing  the  number  of  time  steps  by 
modifying  the  stability  condition  did  not  reduce  the  contour  fluctuations  and  it’s 
likely  the  oscillations  are  an  artifact  associated  with  the  numerical  viscosity 
inherent  in  the  LWAF  scheme.  Recall  TVD  weight  limiting  is  analogous  to 
artificial  viscosity  methods.  Because  the  two-dimensional  code  applies  weight 
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y  (m) 


limiting  at  every  cell  interface  at  every  time  step,  it’s  conceivable  the  artificial 


viscosity  manifests  itself,  as  time  progresses,  as  numerical  error. 


100  x  100  Cells  200  x  200  Cells 


0  0.5  1  1.5  2  0  0.5  1  1.5  2 


x(m)  x  |’m'| 

Figure  5.6.  Cylindrical  Explosion  Density  Profile  Mesh  Refinement 

As  an  additional  verification  of  circular  symmetry,  densities  were 
compared  along  a  line  extending,  radially,  one  meter  outward  from  (1,1)  at  angles 
of  6  =  0,  7r  /  4,  7t/2,  3ir /  4,  and  tv  ,  as  measured  from  the  x-axis.  Figure  5.7  shows 

the  density  profiles  along  each  of  these  lines  and  indicates  nearly  identical  values 
for  each  angle. 

Based  on  these  analyses,  I  conclude  two-dimensional  flows  are  properly 
modeled  using  dimension  splitting  techniques. 
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Figure  5.7.  Mild  Cylindrical  Explosion  Density  Comparison  at  t  =  0.25  sec 
5.3.  Implementation  in  Three  Dimensions 
5.3.1.  Initial  Boundary  Value  Problem 


The  three  dimensional  initial  boundary  value  problem  for  non-linear 


systems  of  hyperbolic  conservation  laws  is 

PDFs  :  U*  +  F(U)X  +  G(U)„  +  H(U)Z  =  0, 
ICs  :  U(x,y,z,0)  =  Xj(°\x,y,z). 


Here  is  a  piecewise  continuous  distribution  of  initial  data  defined  over 

the  spatial  domain  [0,  Lx]  x  [0,  Ly]  x  [0,  LJ  at  t  =  0.  Lastly,  in  three  dimensions, 
boundary  conditions  are  described  along  a  plane  at  each  of  the  six  boundaries 
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using  symmetric  and  transmissive  boundary  conditions,  as  described  in  the  one¬ 


dimensional  case. 

5.3.2.  Discretization  of  the  Euler  Conservation  Equations 

Discretization  of  the  spatial  and  temporal  computational  domain  [0,  Lx]  x 
[o,  L,]  x  [0,  LJ  x  [0,  T]  requires  extension  of  the  two-dimensional  mesh  in  the  in¬ 
direction.  The  spatial  mesh  is  partitioned  into  /  x  J  xK  computational  cells  of 
uniform  volume  Ax  x  Ay  x  Az  such  that 

Ax  = 

I 

Av=Iji  (5.22) 

Az  =  —. 

K 

Each  computational  cell  is  bounded  by  six  faces  i-  1/2,  i  +  1/2,  j  -  1/2,  j  +  1/2, 
k  -  1/2,  k  +  1/2  in  the  x,  y,  and  z  directions,  respectively,  where 

xi- 1/2  =  (i-l)Ax, 

xi+l/2  =  iAxi 

Vj -i/2  =  U  ~  i)Ay, 

7  (5.23) 

2/J+1/2  =  jAv, 

zk- 1/2  =  (fc-l)^, 

*fc+l/2  = 

The  center  of  each  computational  cell  ( i,  j,  k)  is  determined  as 

(Xi,yj,zk)  =  ((*  -  |)Ar,  (j  -  (fc  - 1)^).  (5.24) 
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As  in  the  one  and  two-dimensional  cases,  the  temporal  domain  is  separated  into  n 


non-uniform,  adaptively  determined  time  steps  where 

CcflAy  c_cfiAz 


At  =  min 


C^Ax 


(SfW  (^Lax  m 


(5.25) 


(S'”  )max  ,  (S'-  )max  ,  and  (S%  )max  are  approximations  for  maximum  wave  speeds  in 
the  x,  y,  and  z  directions,  respectively,  during  time  level  n,  and  are  based  upon 


the  material  and  sound  speeds  such  that 

)max  max  j 

In 

;  |  +  ai,  j,k  } 

(Sy  )max  max  | 

fe* 

:|  +  ai,j,k  } 

)max  max  | 

ki, 

k  |  A  az,j,k  j 

(5.26) 


Finally,  as  in  the  one  and  two-dimensional  models,  the  conditions  set  forth 
ultimately  predict  stable  solutions  across  the  entire  computational  domain  during 
each  time  step. 


5.3.3.  Dimension  Splitting 

The  splitting  techniques  executed  in  three  dimensions  are  analogous  to  the 
two-dimensional  procedures.  Equation  (5.21)  is  resolved  as  three  IBVPs  (Toro, 


1999:543) 


PDEs  : 
ICs  : 


U,  +  F(U)j.  =  0. 


u 


n 

i,j,k 


At 


n+1/3 
i,j,k  ’ 


(5.27) 
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jjn+2/ 3 


(5.28) 


PDEs  :  XJt  +  G(U)y  =  0, 


ICs 


tj71+1/  3 


and 


PDEs:  Ut+H(U)z 


ICs  : 


U 


rc+2/3 

i,j,k 


0, 


At. 


•u 


n+1 

i,j,k 


(5.29) 


which  may  be  summarized,  in  operator  notation,  as 

XJl+l  =  Z{At)Y{At)X{At\\jXk).  (5.30) 

As  in  the  two-dimensional  case,  analysis  of  this  method  indicates  first-order 
accuracy  in  time.  However,  second-order  spatial  accuracy  is  achieved,  every 
other  time  step,  by  the  method 

Ugf  =  x^At)Y^At) Z^At) Z^At)Y^At) X^At\ VPk) .  (5.31) 

In  practice,  (5.31)  is  implemented  over  the  interval  [tn,tn+2]  based  upon  the 
following  method.  Given  initial  conditions  ,  boundary  conditions  are 

established  as  symmetric  or  transmissive  and  a  stable  time  step  is  approximated 
via  (5.25).  The  first  of  two  time  steps  begins  with  X^At^  operating  on  XJXk . 
Specifically,  in  the  ^-direction,  .J  x  K  one-dimensional  Riemann  problems  are 
solved,  fluxes  are  calculated  at  each  computational  boundary,  and  temporary 
variable  solutions  are  updated  in  each  cell  as  Uy,*  .  Next,  boundary  conditions 
are  updated  along  the  planes  y  =  0  and  y  =  Ly  and  operates  on  Uijjt  by 
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solving  I  xK  one-dimensional  problems  in  the  ^-direction,  as  described  above, 

^71  +  1 

producing  updated  cell  solutions  Ui,j,k  .  Before  solving  problems  in  the  z- 
direction,  boundary  conditions  are  applied  along  the  planes  z  =  0  and  z  =  Lz  . 

At  last,  operates  on  IXyi;  by  solving  JxJ  one-dimensional  problems  in  the 
^-direction.  At  the  conclusion  of  this  first  time  step,  solutions  are  analyzed 
against  the  CFL  stability  criteria.  If  any  CFL  criterion  is  violated,  then:  the 
interim  solutions  are  discarded;  the  CFL  coefficient  and  time  step  are  re¬ 
calculated;  and  the  time  step  is  reset.  However,  provided  stable  solutions  prevail, 
the  second  time  step  is  executed  over  solution  U^jJ; ,  with  time  step  operator 
ordering  reversed  to  eliminate  directional  bias.  As  in  the  two-dimensional  case, 
subject  to  stability  verification  of  ,  simulation  time  is  advanced  by  2  At  and 

the  second-order  accurate  scheme  is  repeated  until  the  desired  solution  time  is 
achieved. 

5.3.4 •  Special  Considerations 

Three-dimensional  splitting  facilitates  solving  the  Euler  conservation 

equations,  in  the  x,  y,  and  z  directions,  as 

if  pu 
P 

pu  pu 2  +  p 

E  +  E(u  +  p)  =  0,  (5.32) 

pv  pvu 

pw  pwu 

t  1  IX 
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and 


p 

pv 

pv 

e 

to 

+ 

E 

+ 

E(v  +  p ) 

pu 

puv 

pw 

t 

pwv 

p 

pw 

pw 

pw 2  +  p 

E 

+ 

E(w  +  p) 

pu 

puw 

pv 

t 

pvw 

=  0, 


=  0. 


(5.33) 


(5.34) 


As  in  the  two-dimensional  case,  tangential  velocities,  their  associated  momentum 
fluxes,  and  the  energy  flux  require  special  consideration.  Within  system  (5.32),  u 
represents  the  normal  component  of  material  velocity,  while  v  and  w  represents 
the  tangential,  or  shear,  components.  The  normal  velocity  component  in  system 
(5.33)  is  v  while  u  and  w  signify  tangential  components.  Finally,  in  (5.34),  w 
denotes  the  normal  velocity  component,  while  u  and  v  correspond  to  tangential 
components. 

Solutions  to  the  aforementioned  systems  of  equations  that  apply  TVD  flux 
limiter  functions  produce  the  dissipative  wave  structure  of 

Figure  5.8.  As  in  the  two-dimensional  case,  normal  and  tangential  flux 
calculations  are  determined  separately.  Similarly,  weighted  average  states  of 
density,  normal  velocity,  and 
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Figure  5.8.  Three-Dimensional  TVD  Flux  Limiter  Wave  Structure 

pressure  are  determined  using  equations  (4.43)  and  (4.44),  and  their 
corresponding  fluxes  are  subsequently  computed.  Calculation  of  the  tangential 
velocity  weighted  average  states  and  their  resultant  tangential  momentum  fluxes 
is  also  analogous  to  the  two-dimensional  case.  Specifically,  if  Criand  VT2  denote 
the  two  tangential  velocity  components,  then  the  weighted  average  states  of  these 
elements  are  calculated  as 


Cry+l/2  —  ^5,LWAF^T1,1+  ^6,LWAF^T1,4 

VT2,i+l/2  =  Pt,LWAfVt2.\  +  f^8,LWAF^T2,A 


(5.35) 


where  P^lwaf  anfi  WAF are  given  by  equations  (5.17)  and  (5.18)  and 

P 7,lwaf  =  (1  +  (f)5Sign(c5)) , 

/%  ,LWAF  =  (1  —  4>5^^9n(c5  )  )  • 
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(5.36) 


Conveniently,  the  eigenvalue  structure  of  solutions  to  (5.32)  -  (5.34)  dictates 
c5  =  c4  =  °2  an(l  -  like  04  ,  the  weight  limiter  function  05  is  calculated  from  its 
tangential  velocity  gradients  AVT2as  in  equation  (5.19).  Hence,  calculation  of  all 
components  required  to  solve  the  three-dimensional  time  dependent  Euler 
conservation  equations  is  complete. 

5.3.5.  Sod  Test  in  Three  Dimensions 

The  benchmark  shock  tube  test  is  re-visited  here  for  the  final  time,  to 

verify  correct  modeling  of  planar  flows  in  three  dimensions.  The  three 

dimensional  Sod  test  was  executed  over  the  cubic  domain  [0,  lm]  x  [0,  lm]  x  [0, 

lm]  consisting  of  100  cells  in  the  x,  y,  and  z  directions  and  a  diaphragm, 

positioned  along  the  plane  x  =  0.5  m.  Initial  conditions  were  established  as 

WL  =  (1.000  kg/m3, 0.0 m/sec,  O.Om/sec,  0.0m/sec,1.0N/m2)  x  <  0.5m, 

W(x,  y,  z,  0)  =  •  (5.37) 

WR  =  (0.125 kg/m3, O.Om/sec,  O.Om/sec,  O.Om/sec, 0.1N/m2)  x  >  0.5m 

and  solutions  developed  over  the  time  interval  t  =  [0,  0.25  sec].  The  density 
profile  at  t  =  0.25  sec,  through  the  plane  z  =  0.5m  is  depicted  in  Figure  5.9  and 
is  consistent  with  expected  results.  As  in  the  two-dimensional  case,  correct 
modeling  of  planar  flows  was  verified  by  comparing  data  extracted  from  the  first 
set  of  cells  in  the  ^-direction  against  the  exact  solution  and  the  second-order 
accurate  solution  obtained  by  the  LWAF  method  in  a  single  dimension.  Again, 
the  second-order  accurate  solutions  appeared  graphically  identical  for  both  the 
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one  and  three-dimensional  cases,  and  I  conclude  that  the  developed 
computational  model  accurately  models  planar  flows. 


Figure  5.9.  Three-Dimensional  Sod  Density  Profile  (z  =  0.5m)  at  t  =  0.25  sec 
5.3.6.  Mild  Spherical  Shock  Test 

To  verify  that  the  developed  code  correctly  models  flows  in  three  dimensions,  the 
Euler  equations  were  solved  over  the  cubic  domain  [0,  2m]  x  [0,  2m/  x  [0,  2m] 
discretized  as  a  computational  mesh  of  100  cells  in  the  x,  y,  and  z  directions. 
Initial  conditions  consisted  of  a  spherical  region  with  radius  0.4m,  centered  at 
(1,1,1),  with  primitive  variables  of  WL and  TE^from  (5.37)  representing  states 
inside  and  outside  the  sphere,  respectively  (Toro,  1999:587).  Figure  5.10  shows 


90 


the  density  profile  through  the  plane  z  =  1  m  at  t  =  0.25  sec  and  indicates,  as 
expected,  a  solution  which  exhibits  a  spherical  shock,  contact  surface,  and 
rarefaction. 


Figure  5.10.  Spherical  Explosion  Density  Profile  (z  =  lm)  at  t  =  0.25  sec 

As  in  the  two-dimensional  case,  the  density  profile  indicates  the  presence 
of  potential  numerical  artifacts  in  the  region  between  the  rarefaction  and  contact 
discontinuity.  Accordingly,  a  similar  analysis  of  contour  plots  and  mesh 
refinement  was  conducted.  Increasing  cells  in  all  three  directions  reduced  the 
oscillations  in  the  region  of  concern  and  yielded  improved  circular  symmetry. 
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Furthermore,  densities  compared  along  a  line  extending,  radially,  one  meter 
outward  from  (1,1,1),  in  the  plane  z  =  1  m,  at  angles  of 

6  =  0,  7r  /  4,  7t/2,  3ir  /  4,  and  7r  ,  as  measured  from  the  rr-axis.  Figure  5.11  shows 

the  density  profiles  along  each  of  these  lines  and  indicates  nearly  identical  values 
for  each  angle. 
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Figure  5.11.  Spherical  Explosion  Density  Comparison  (z  =  1)  at  t  =  0.25  sec 

To  confirm  spherical  symmetry,  similar  analyses  were  conducted  through 
planes  y  =  1  m  and  x  =  1  m.  Figure  5.12  and  Figure  5.13  show  the  density 
profiles  through  the  planes  y  =  1  m  and  x  =  1  m,  respectively,  at  t  =  0.25  sec 
and  are  consistent  with  the  solution  presented  in  Figure  5.10.  Densities  were 
compared  along  a  line  extending,  radially,  one  meter  outward  from  (1,1,1),  in  the 
plane  y  =  1  m,  at  angles  0  =  0,  7r  /  4,  7t/2,  37t  /  4,  and  7r ,  as  measured  from  the  x- 

axis.  Figure  5.14  shows  these  density  profiles  along  each  of  these  lines  and 
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indicates  nearly  identical  values  for  each  angle.  Finally,  Figure  5.15  depicts 
identical  densities  plotted  along  a  one  meter  radial  line,  extended  from  (1,1,1),  in 
the  plane  x  =  1  m,  at  angles  17  =  0,  ir  /  4,  ir  /  2,  37r  /  4,  and  tv  ,  measured  from  the 

y- axis.  Based  on  these  analyses,  I  conclude  multi-dimensional  flows  are  properly 
modeled  using  dimension  splitting  techniques. 
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Figure  5.12.  Spherical  Explosion  Density  Profile  ( y  =  lm)  at  t  =  0.25  sec 
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Figure  5.13.  Spherical  Explosion  Density  Profile  ( x  =  lm)  at  t  =  0.25  sec 
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Figure  5.14.  Spherical  Explosion  Density  Comparison  ( y  —  lm)  at  t  —  0.25  sec 
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Figure  5.15.  Spherical  Explosion  Density  Comparison  ( x  =  lm)  at  t  =  0.25  sec 


6.  Conclusion 


The  principal  endeavor  of  this  research  was  development  of  a  three- 
dimensional  hydrodynamic  shock  code  to  solve  the  Euler  conservation  equations, 
while  the  secondary  objective  included  applying  the  developed  code  to  model  air 
blast  propagation. 

To  that  end,  a  one-dimensional,  second-order  accurate  computational 
model  was  methodically  developed,  tested,  and  validated.  Godunov’s  scheme 
executed  with  exact  and  approximate  Riemann  solvers,  and  applied  to  varying 
shock  tube  tests,  accurately  predicted  wave  structure  and  locations.  Although 
both  solvers  yielded  comparable  accuracy,  the  adaptive,  approximate  solver 
improved  computational  efficiency  by  nearly  36%.  Application  of  E.F.  Toro’s 
weighted  average  flux  method  produced  second-order  accuracy,  but  introduced 
false  fluctuations  near  discontinuous  regions.  These  oscillations  were  effectively 
reduced  via  TVD  weight  limiting  techniques,  where  van  Leer’s  limiter  most 
favorably  controlled  behavior  in  the  presence  of  shocks.  The  fully  developed 
second-order  accurate  model  was  tested  and  verified  against  the  ARL  57  cm 
shock  tube  test,  and  agreed  with  the  published  experimental  data  within  5.9%. 

Following  validation  in  a  single  dimension,  the  developed  code  was 
extended  into  two  and  three  dimensions.  During  both  implementations,  output 
from  the  multi-dimensional  models  was  compared  against  Sod’s  shock  test,  to 
verify  accurate  modeling  of  planar  flows,  and  evaluated  against  mild  cylindrical 
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and  spherical  shock  tests,  presented  by  Toro,  to  confirm  correct  flow 
representation  in  two  and  three  dimensions. 

The  end  result  of  this  research  effort  is  a  fully  developed  three-dimensional 
air  blast  propagation  model  that  may  be  implemented,  using  parallel 
computations,  to  permit  a  more  accurate  and  efficient  exploration  of  air  blast 
propagation. 
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diminishing  (TVD)  flux  and  weight  limiting  schemes  are  added  to  reduce  fluctuation  severity.  Finally,  the  fully  developed  one¬ 
dimensional  code  is  validated  against  experimental  data,  and  extended  into  two  and  three  dimensions  via  dimension-splitting 
techniques. 
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